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Abstract 

In  this  thesis  the  analysis  of  natural  ice  events  is  carried  out  based  on  direct  measure¬ 
ments  of  ice-borne  seismo-acoustic  waves  generated  by  ice  fracturing  processes.  A 
major  reason  for  studying  this  phenomenon  is  that  this  acoustic  emission  is  a  signifi¬ 
cant  contributor  to  Arctic  ocean  ambient  noise.  Also  the  Arctic  contains  rich  mineral 
and  oil  resources  and  in  order  to  design  mining  facilities  able  to  withstand  the  harsh 
environmental  conditions,  we  need  to  have  a  better  understanding  of  the  processes 
of  sea  ice  mechanics.  The  data  analyzed  in  this  thesis  were  collected  during  the  Sea 
Ice  Mechanics  Initiative  (SIMr94)  experiment  which  was  carried  out  in  the  spring  of 
1994  in  the  Central  Arctic. 

One  of  the  contributions  of  this  thesis  was  the  determination  of  the  polarization 
characteristics  of  elastic  waves  using  multicomponent  geophone  data.  Polarization 
methods  are  well  known  in  seismology,  but  they  have  never  been  used  for  ice  event 
data  processing.  In  this  work  one  of  the  polarization  methods  (so  called  Motion  Prod¬ 
uct  Detector  method)  has  been  successfully  applied  for  localization  of  ice  events  and 
determination  of  polarization  characteristics  of  elastic  waves  generated  by  fractur¬ 
ing  events.  This  application  demonstrates  the  feasibility  of  the  polarization  method 
for  ice  event  data  processing  because  it  allows  one  to  identify  areas  of  high  stress 
concentration  and  “hot  spots”  in  ridge  building  process. 

The  identification  of  source  mechanisms  is  based  on  the  radiation  patterns  of  the 
events.  This  identification  was  carried  out  through  the  analysis  of  the  seismo-acoustic 
emission  of  natural  ice  events  in  the  ice  sheet.  Previous  work  on  natural  ice  event 
identification  was  done  indirectly  by  analyzing  the  acoustic  energy  radiated  into  the 
water  through  coupling  from  elastic  energy  in  the  ice  sheet. 

After  identification  of  the  events,  the  estimation  of  the  parameters  of  fault  pro¬ 
cesses  in  Arctic  ice  is  carried  out.  Stress  drop,  seismic  moment  and  the  type  of  ice 
fracture  are  determined  using  direct  near-field  measurements  of  seismo-acoustic  sig¬ 
nals  generated  by  ice  events.  Estimated  values  of  fracture  parameters  were  in  good 


agreement  with  previous  work  for  marginal  ice  zone. 

During  data  processing  the  new  phenomenon  was  discovered:  “edge  waves” ,  which 
are  waves  propagating  back  and  forth  along  a  newly  opened  ice  lead.  These  waves 
exhibit  a  quasi-periodic  behavior  suggesting  some  kind  of  stick-slip  generation  mecha¬ 
nism  somewhere  along  the  length  of  the  lead.  The  propagation  characteristics  of  these 
waves  were  determined  using  seismic  wavenumber  estimation  techniques.  In  the  low- 
frequency  limit  the  dispersion  can  be  modeled  approximately  by  an  interaction  at  the 
lead  edges  of  the  lowest  order,  antisymmetric  modes  of  the  infinite  plate. 

Thesis  Supervisor:  Henrik  Schmidt 

Title:  Professor,  Massachusetts  Institute  of  Technology 
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Chapter  1 


Introduction 


1.1  Motivation 

For  several  years  there  has  been  a  great  interest  in  arctic  ice  mechanics  research 
even  after  the  arctic  region  lost  its  strategic  importance  for  military  applications. 
One  of  the  reasons  is  that  the  ice  cover  can  serve  as  a  sensitive  indicator  of  changes 
in  global  climate  (e.g.  the  greenhouse  effect).  Another  reason  is  that  the  Arctic 
contains  rich  mineral  and  oil  resources  and  in  order  for  mining  facilities  to  be  able  to 
withstand  harsh  environmental  conditions,  we  need  to  have  a  better  understanding 
of  the  processes  of  the  sea  ice  mechanics.  This  knowledge  can  also  help  in  improving 
navigation  conditions  for  ships  in  arctic  seas. 

1.2  Objectives 

An  example  of  this  interest  is  the  Sea  Ice  Mechanics  Initiative  (SIMI)  experiment 
carried  out  in  a  joint  effort  by  several  institutions  in  Spring  of  1994.  The  objective  of 
SIMI  was  to  obtain  basic  physical  understanding  of  the  highly  non-linear  process  by 
which  environmental  forcing  (including  atmospheric  conditions,  oceanic  currents,  heat 
flux  etc.)  leads  to  formation  of  the  dramatic  features  characterizing  the  macroscopic 
Arctic  ice,  according  to  ([52], [45]).  This  objective  required  a  method  of  monitoring  the 
mechanical  processes  in  ice  that  could  determine  the  development  of  stress  and  strain 
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distributions  in  space  as  well  as  time,  providing  both  coverage  and  resolution.  This 
requirement  suggested  acoustic  and  seismo-acoustic  remote  sensing  as  an  attractive 
and  feasible  option.  One  such  remote  sensing  method  is  based  on  the  detection  and 
analysis  of  the  stress  waves  generated  in  the  ice  and  the  water  column  by  fracture 
formation  in  the  ice  cover.  When  a  crack  is  formed,  the  stresses  are  released  over 
the  crack  surface,  producing  stress  waves  radiated  in  a  pattern  characteristic  for 
that  particular  fracture  process.  The  parameters  describing  the  fracture  process  can 
therefore  be  determined  by  inversion  of  the  radiated  stress  wave  field. 

The  traditional  experimental  setup  for  remote  sensing  of  ice  events  used  the  hor¬ 
izontal  hydrophone  array  to  monitor  the  fracture  development  in  the  far-field.  This 
approach  does  not  provide  sufficient  spatial  resolution  to  invert  for  all  fracture  pa¬ 
rameters  describing  the  ice  mechanical  development.  To  achieve  this  the  sensors  need 
to  be  deployed  in  the  vicinity,  or  near- field,  of  the  ice  fracture  events. 

Recent  experiments  ([33])  have  demonstrated  a  resolution  superiority  of  three- 
component  geophones  compared  to  hydrophones  for  near-field  observation  of  seismo- 
acoustic  propagation  in  the  ice  sheet.  Additionally,  the  geophones  allow  direct  mea¬ 
surement  of  the  characteristics  of  the  wave  propagation  in  the  ice,  whereas  hy¬ 
drophones  give  information  only  about  the  part  of  the  event  energy  which  was  trans¬ 
ferred  to  acoustic  waves  in  the  water  due  to  the  interaction  at  the  ice-water  interface. 
As  a  result,  the  ice  crack  characteristics  appear  in  hydrophone  data  only  as  second 
order  effects.  Also,  the  vector  nature  of  the  data  recorded  by  three-component  geo¬ 
phones  provides  substantially  richer  information  than  hydrophones  about  the  char¬ 
acteristics  of  the  waves  in  the  ice  sheet. 

In  the  following  I  will  analyze  and  model  the  geophone  data  collected  by  the 
joint  MIT/WHOI  group  during  SIMI-94  experiment  in  order  to  address  some  of  the 
objectives  of  Sea  Ice  Mechanics  Initiative.  The  specific  objectives  of  this  analysis  will 
be  identification  of  ice  events  and  estimation  of  their  important  parameters  from  the 
elastic  waves  generated  by  the  events  in  the  ice  cover. 
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1.3  Previous  work 


This  review  of  the  previous  work  is  organized  along  the  following  lines:  first,  I  will 
discuss  two  directions  in  the  study  of  the  Arctic  ambient  noise,  then  I  will  review 
literature  devoted  to  the  application  of  earthquake  seismology  methods  to  the  analysis 
of  ice  event  data. 

The  study  of  the  Arctic  ambient  noise  usually  follows  two  distinct  directions.  One 
direction  is  to  correlate  different  characteristics  of  ambient  noise  to  environmental 
parameters  such  as  heat  fiux,  temperature  and  wind.  The  other  direction  is  to  study 
individual  ice  events  which  form  the  average  ambient  noise  in  Arctic.  In  that  approach 
the  implicit  assumption  is  that  the  ambient  noise  in  Arctic  is  dominated  by  the 
mechanical  processes  in  the  ice  cover  [8]. 

The  first  direction  started  with  work  by  Milne  and  Ganton  [35].  They  identified 
the  major  source  of  ambient  noise  for  shore-fast  ice  during  Spring  and  Winter  as 
thermal  tensile  cracks,  while  during  Summer  -  as  relative  motion  of  ice  floes.  In 
their  next  work  [12]  for  Winter  season  they  identified  two  different  spectral  regions 
for  ambient  noise:  lower  frequencies  (between  200  and  800  Hz)  showed  an  impulsive 
character  of  radiation  with  noise  level  declining  rapidly  with  increase  of  temperature; 
the  other  frequency  range  (between  1  kHz  and  10  kHz)  had  Gaussian  character  and 
was  attributed  to  wind  forcing  due  to  dependence  on  wind  speed. 

Makris  and  Dyer  [29]  showed  that  long-term  variations  of  low-frequency  ambient 
noise  under  pack  ice  of  the  central  Arctic  Ocean  correlated  highly  with  composite 
measures  of  stress  applied  to  the  ice  by  wind,  current,  and  drift.  These  compos¬ 
ites  were  identified  as  the  horizontal  ice  stress  and  the  stress  moment,  and  were 
derived  from  meteorological  and  oceanographic  data  observed  simultaneously  with 
the  noise.  They  also  concluded  that  atmospheric  cooling  (a  known  high  correlate  of 
mid-frequency  noise  under  the  ice)  was  not  important  at  low  frequencies.  In  contrast 
to  the  results  for  the  central  Arctic,  they  showed  [30]  that  in  the  marginal  ice  zone 
ice  stress,  ice  moment,  and  wind  stress  magnitude  were  poor  correlates  for  the  low- 
frequency  ambient  noise.  They  concluded  that  for  that  Arctic  region  the  primary 
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correlate  of  the  noise  was  surface  gravity  wave  forcing,  with  the  most  likely  mech¬ 
anisms  for  sound  generation  being  the  flexural  floe  failure  and  unloading  motion, 
within  a  few  kilometers  of  the  ice  edge. 

The  study  of  event  source  mechanisms  was  started  by  Milne.  He  considered  ther¬ 
mal  tensile  cracks  as  the  noise  source  mechanism  [34]  and  found  that  the  ambient  noise 
level  in  the  150-300  Hz  band  correlated  well  with  the  changes  in  air  temperature.  In 
[36]  he  suggested  that  the  dominating  mechanisms  of  ambient  noise  generation  may 
vary  with  the  Arctic  region.  This  assumption  was  later  confirmed  in  other  exper¬ 
iments.  Lewis  et  al.  [26], [27]  corrected  some  of  the  conclusions  of  the  theory  for 
thermal  cracking  in  ice  by  Milne  et  al.  They  showed  that  heat  flux,  rather  than  air 
temperature,  was  the  main  factor  in  ice  fracturing  under  thermal  tension  [26].  Their 
other  main  contribution  to  the  study  of  event  mechanisms  was  the  calculation  of  the 
ice  stress  due  to  the  thermal  effect  [27]. 

Pritchard  [40]  suggested  that  the  energy  dissipated  during  the  ridging  process  was 
the  proper  measure  of  the  ambient  noise  source  level.  Numerical  simulations  showed 
that  a  significant  amount  of  energy  and  ambient  noise  could  be  explained  that  way. 
The  model  was  also  physically  attractive  and  properly  explained  lack  of  noise  when 
winds  were  high  but  the  ice  was  strong  enough  to  resist  ridging. 

The  theoretical  modeling  of  ice  cracking  events  as  a  source  of  acoustic  energy 
radiated  into  the  water  was  started  by  Stein  [50].  He  used  a  monopole  source  in 
the  ice  plate  to  model  the  tensile  ice  cracking.  His  was  an  ad  hoc  radiation  model 
appropriate  only  for  events  whose  openings  were  small  compared  to  ice  thickness.  For 
more  complex  events  a  monopole  source  should  be  combined  with  a  quadrupole  source 
as  suggested  by  Langley  [23],  [24],  [25].  Kim  [22]  used  a  computational  approach  to 
the  crack  radiation  problem.  He  developed  an  analytical  and  numerical  model  of  the 
elastic  wave  field  in  range  independent  elastic  environments  for  various  seismic  source 
mechanisms  including  shear  and  tensile  cracks.  Using  compact  source  representation 
with  fault  surface  in  arbitrary  direction,  he  derived  the  basic  Green’s  function  solution 
for  propagation  in  stratified  elastic  media  which  was  incorporated  into  a  numerical 
model.  Kim  extended  that  solution  to  include  more  complete  cracking  mechanisms. 
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like  non-compact  and  moving  cracks.  He  also  studied  some  effects  of  anisotropy  on 
the  acoustic  emission. 

Dyer  [7],  [8]  discussed  several  possible  event  mechanisms  of  noise  generation  in 
the  low-  and  mid-frequency  ranges  of  the  noise  spectrum  with  the  ridge  unloading 
model  among  them.  He  also  introduced  a  slip  model  for  ice  fracture.  Chen  [4]  used 
these  models  in  her  study  of  ice  event  mechanisms  in  the  Marginal  Ice  Zone.  The 
application  of  similar  ideas  to  the  central  Arctic  region  was  carried  out  by  Stamoulis 
in  her  thesis  [49].  She  also  introduced  event  identification  by  their  radiation  patterns 
into  the  study  of  events  mechanisms. 

The  above  studies  were  mostly  concerned  with  the  sound  radiated  into  the  water 
by  ice  events.  The  more  direct  way  to  study  ice  events  mechanisms  is  to  investigate 
the  elastic  field  produced  by  the  event  in  the  ice  sheet  itself.  This  could  be  achieved  by 
using  geophones,  sensors  that  are  coupled  directly  to  the  ice  by  freezing,  and  capable 
of  measuring  the  particle  velocity  of  ice  motion  directly.  One  of  the  first  studies  of 
ice  events  in  the  Arctic  using  geophones  was  conducted  by  Hunkins  in  1960  [18].  He 
successfully  managed  to  isolate  the  flexural  waves  in  ice  by  investigating  the  sense  of 
particle  orbital  motion  in  the  wave  near  the  ice  surface.  For  this  purpose  he  plotted 
the  output  of  the  vertical  axis  of  the  sensor  vs  the  output  of  one  of  the  horizontal  axes. 
His  diagrams  could  not  show  the  true  particle  motion  due  to  frequency  dependence 
of  the  phase  shift.  His  results  showed  the  elongated  horizontal  ellipse  with  retrograde 
direction  of  motion  just  like  it  was  predicted  by  Sato  [42], [43].  According  to  Hunkins, 
only  flexural  waves  gave  the  clear  picture  of  the  particle  orbital  motion,  other  waves 
just  produced  a  jumble  on  the  hodographs,  whereas  later  [33]  longitudinal  waves 
have  been  also  identified.  In  some  situations  this  method  could  help  to  distinguish 
different  wave  types  by  their  planes  of  polarization.  Note  that  the  above  hodograph 
method  could  only  be  applied  to  data  from  multi-axis  sensors  such  as  three-component 
geophones. 

Relatively  recently,  Yang  and  Giellis  used  an  array  of  geophones  for  the  determi¬ 
nation  of  dispersion  characteristics  of  flexural  wave  in  ice  generated  by  striking  sledge 
hammers  against  a  steel  pipe  frozen  into  the  ice  during  their  Arctic  experiment  in 


22 


1988  [54],  [55].  However,  they  used  only  the  vertical  component  of  velocity  recorded 
by  the  geophones  in  the  data  processing.  The  data  processing  scheme  made  use  of 
the  Choi- Williams  display  [5]  for  selected  vertical-axis  geophones  in  which  the  signal 
intensity  was  plotted  as  a  function  of  time  and  frequency.  Then  those  displays  were 
converted  into  group-velocity  versus  frequency  plots  using  the  known  receiver  range 
and  measured  travel  time.  Finally,  converted  plots  (associated  with  nine  vertical-axis 
geophones  at  different  ranges)  were  summed  in  intensity  (incoherent  processing). 

B.E.  Miller  ([33],  [32])  used  geophone  data  to  precisely  locate  ice  events  in  the 
PRUDEX-87  experiment.  He  found  that  the  geophones  give  much  more  precise  es¬ 
timates  of  locations  for  experimental  shots  than  hydrophones  do.  He  also  used  that 
data  to  estimate  the  location  of  an  ice  ridge.  In  the  geophone  time  series  generated 
by  under-ice  detonations,  he  found  not  only  the  expected  longitudinal  and  flexural 
waves  in  the  ice  plate,  but  also  unexpected  horizontally  polarized  transverse  (SH) 
waves.  The  need  to  determine  the  origin  of  the  SH  wave  from  the  received  time  series 
additionally  highlighted  the  dramatic  superiority  of  geophones  over  hydrophones  in 
this  application.  Miller  also  conducted  the  inversion  of  the  geophone  data  for  the 
ice  sheet’s  low-frequency  elastic  parameters  initially  by  modeling  the  ice  as  a  single 
homogeneous  isotropic  plate  using  SAFARI.  After  that  a  modified  stationary  phase 
approach  was  used  to  extend  SAFARI  modeling  to  invert  for  the  parameters  of  two 
ice  half-plates  simultaneously. 

The  stress  relief  in  arctic  ice  have  features  similar  to  that  of  the  stress  relief  in  the 
earth’s  crust.  Hence,  earthquake  mechanics  should  provide  a  useful  starting  point  for 
examining  seismo-acoustic  signals  radiated  by  the  ice  crack  [10].  By  far,  the  most 
common  types  of  cracks  in  ice  mechanics  are  tensile  cracks,  dip-slip  and  strike-slip 
(both  with  dip  angle  5  0“  [22]).  For  these  kinds  of  cracks  the  general  model  was 

adapted  by  D.  Farmer  and  Y.  Xie  in  [10].  By  analogy  with  results  of  theoretical 
models  of  propagating  earthquake  cracks,  they  assumed  that  the  source  moves  with 
the  developing  crack.  The  authors  adapted  a  concept  introduced  in  seismology  by 
Haskell  and  modeled  the  source  as  a  sinusoidally  roughened  ramp  function,  thereby 
inducing  a  significant  fraction  of  the  energy  in  a  primary  or  base-band  lobe  (200-300 
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Hz).  The  frequency  of  this  lobe  is  governed  by  the  propagation  speed  and  by  the 
coherent  length  of  a  crack  segment.  The  high-frequency  peak  (pi5  kHz)  in  this  model 
is  related  to  the  vertical  length  scale  of  the  crack.  In  this  approach  the  calculation 
of  event  parameters  requires  an  assumption  regarding  the  rupture  velocity  Vr-  Using 
[31]  one  gets 


Vr{max)  ^  0.63/3,  (1.1) 

where  is  the  shear  wave  speed.  After  that  the  crack  length  L  is  obtained  from  the 
base-band  frequency  /&  as 


^  Vr 

This  actually  gives  the  length  of  the  coherent  crack  segment  (segment  on  which  crack 
can  be  treated  as  a  fault  propagation  having  constant  speed). 

The  other  similarity  of  seismology  to  seismo-acoustics  is  the  use  of  similar  sen¬ 
sors.  Both  fields  use  devices  that  measure  particle  velocity  of  the  medium  quite 
extensively,  though,  unlike  seismo-acoustics,  in  seismology  an  array  of  other  instru¬ 
ments  is  also  deployed  (among  them  accelerometers  and  devices  that  measure  particle 
displacement  directly).  However,  until  recently  (about  mid-80s)  seismologists  did  not 
really  use  arrays  of  three-component  sensors.  They  limited  their  attention  to  arrays 
of  one-component  sensors  or  single  three-component  sensors.  To  my  knowledge,  the 
paper  [19]  by  A.  Jurkevics  pioneered  the  use  of  the  array  of  three-component  sensors 
in  seismology.  The  analysis  technique  in  this  research  was  based  on  a  time-domain 
algorithm  originally  proposed  by  Flinn  in  1965.  In  this  algorithm  the  polarization 
properties  of  the  seismic  waves  were  computed  within  sliding  time  windows  by  solv¬ 
ing  the  eigenproblem  for  the  covariance  matrix  for  three  components  of  motion  of 
each  geophone  separately.  This  technique  was  extended  to  multiple  three-component 
sensors  in  an  array  configuration  by  averaging  covariance  matrices  for  the  different 
sensors.  In  this  case  a  1/M  reduction  in  the  estimation  variance  was  obtained  (M  is 
the  number  of  sensors),  when  the  noise  and  local  scattering  effects  are  uncorrelated. 
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1.4  Thesis  contributions 


•  Determination  of  polarization  characteristics  of  the  elastic  waves  using  multi- 
component  geophone  data. 

Polarization  methods  are  well  known  in  seismology,  but  they  never  have  been 
used  for  ice  event  data  processing.  In  this  work  one  of  the  polarization  methods 
(so  called  Motion  Product  Detector  method)  has  been  successfully  applied  for 
localization  of  ice  events  and  determination  of  polarization  characteristics  of 
the  elastic  waves  generated  by  these  events.  This  application  demonstrated 
the  feasibility  of  using  the  polarization  method  for  ice  event  data  processing, 
because  it  allows  one  to  identify  areas  of  high  stress  concentration  and  “hot 
spots”  in  ridge  building  process. 

Furthermore,  the  polarization  characteristics  of  elastic  waves  could  be  used  for 
establishing  the  physics  of  the  processes  in  the  ice  sheet  which  lead  to  generation 
of  ice  events. 

•  The  discovery  of  the  new  phenomenon:  “edge  waves” ,  which  are  waves  propa¬ 
gating  back  and  forth  along  a  newly  opened  ice  lead. 

These  waves  exhibit  a  quasi-periodic  behavior  suggesting  some  kind  of  stick-slip 
generation  mechanism  somewhere  along  the  length  of  the  lead.  The  propaga¬ 
tion  characteristics  of  these  waves  were  determined  using  seismic  wavenumber 
estimation  techniques.  In  the  low-frequency  limit  the  dispersion  can  be  modeled 
approximately  by  an  interaction  at  the  lead  edges  of  the  lowest  order,  antisym¬ 
metric  modes  of  the  infinite  plate. 

•  Identification  of  source  mechanisms  for  events  based  on  their  radiation  patterns. 

This  identification  was  carried  out  using  near- field  seismo-acoustic  emission  of 
natural  ice  events.  As  far  as  I  know,  this  has  never  been  done  before  directly. 
Previous  work  on  natural  ice  event  identification  was  done  indirectly:  using 
acoustic  energy  radiated  into  the  water  through  coupling  of  into  water  column 
of  elastic  energy  in  the  ice  sheet  radiated  by  ice  events. 
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•  Estimation  of  parameters  of  fault  processes  in  Arctic  ice. 

Stress  drop,  seismic  moment  and  the  type  of  ice  fracture  were  determined  us¬ 
ing  direct  near-field  measurements  of  seismo-acoustic  signals  generated  by  ice 
events.  Obtained  values  of  fracture  parameters  were  in  good  agreement  with 
those  obtained  by  C.F.  Chen  [4]  for  marginal  ice  zone. 

1.5  Thesis  organization 

First,  in  Chapter  2 1  give  a  short  overview  of  the  SIMI-94  experiment,  data  description 
and  initial  stages  of  data  processing:  preprocessing,  event  detection  and  localization 
using  traditional  least-square  fit  to  interarrival  times.  At  the  end  of  this  chapter  I 
show  some  results  for  flexural  wave  arrival.  In  Chapter  3,  the  polarization  processing 
of  the  geophone  data  is  described.  This  method  allows  one  to  take  advantage  of 
the  multi-component  nature  of  the  geophone  data.  Chapter  4  is  devoted  to  the  new 
phenomenon  found  in  the  geophone  data  collected  in  the  vicinity  of  an  ice  lead,  the 
so-called  “edge  waves”:  waves  propagating  back  and  forth  along  an  ice  lead.  In 
Chapter  5  I  describe  the  results  of  event  source  mechanisms  identification  based  on 
radiation  patterns.  Chapter  6  is  devoted  to  event  parameter  estimation.  The  major 
parameters  estimated  are  the  seismic  moment  and  stress  drop  in  the  fault.  Finally, 
in  Chapter  7,  I  present  a  summary  of  the  thesis  and  suggestions  for  future  work. 
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Chapter  2 


Sea  Ice  Mechanics  Initiative 
(SIMI)  Experiment  overview  and 
data  preparation  procedure 

2.1  Overview  of  the  SIMI  experiment 

As  part  of  the  Office  of  Naval  Research  (ONR)  Sea  Ice  Mechanics  Initiative  (SIMI),  a 
real-time  monitoring  and  processing  program  for  acoustic  emission  from  ice  fracture 
and  ridge-building  events  was  established.  A  wide  aperture  horizontal  hydrophone 
array  was  used  in  combination  with  a  vertical  line  array  to  record  the  acoustic  signals, 
which  were  then  passed  through  a  focused  beamformer  for  real-time  generation  of  ice 
seismicity  maps.  A  number  of  rapidly  deployable  geophone  arrays  were  used  in  active 
zones  to  measure  the  acoustic  emissions  in  the  near  field  for  detailed  seismic  event 
analysis. 

The  primary  objective  of  the  ONR  Sea  Ice  Mechanics  Initiative  was  to  develop 
a  fundamental  understanding  of  the  mechanical  behavior  of  sea  ice  undergoing  envi¬ 
ronmental  forcing  due  to  currents,  wind  and  heat  flux.  The  several  Arctic  acoustics 
experiments  of  the  past  have  made  it  clear  that  ice  fracturing  is  the  dominant  source 
of  ambient  noise,  in  turn  suggesting  that  this  phenomenon  forms  the  major  compo- 
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SIMI-94  Camp  Location 


nent  of  the  physical  processes  leading  from  environmental  forcing  to  the  development 
of  macroscopic  ice  features  such  as  ridges,  leads  and  raftings.  On  this  background,  the 
specific  objective  of  the  MIT/WHOI  SIMI  effort  was  to  develop  a  basic  understanding 
of  these  ice  fracturing  processes  and  their  relation  to  the  environmental  forcing. 

The  main  experiment  of  the  SIMI  effort  was  conducted  in  the  Spring  of  1994 
approximately  400  km  north  of  Prudhoe  Bay,  Alaska  (see  Fig.  2-1).  The  GPS  coor¬ 
dinates  of  the  camp  on  April  22,  1994  were  73°00'56"Ar  and  149°53'55"kF. 

2.2  SIMI  arrays  and  systems 

2.2.1  Surveillance  Hydrophone  Array 

As  a  mechanism  for  monitoring  and  localizing  ice  activity  during  the  Spring  SIMI  field 
experiment,  MIT/WHOI  deployed  two  hydrophone  arrays,  comprising  a  seismicity 
‘surveillance’  capability  approximately  300m  from  the  main  camp  area.  A  vertical 
line  array  (VLA)  of  hydrophones,  having  32  sensors  linearly  spaced  from  a  depth 
62  to  279  meters  was  deployed  at  the  apex  of  the  site  (the  origin  of  the  coordinate 
system  on  the  plot)  shown  in  Fig.  2-2.  Surrounding  this  APEX  was  a  horizontal 
array  consisting  of  a  12  hydrophone  circular  deployment  augmented  by  12  additional 
sensors  on  nominal  north-south  and  east-west  legs  240  m  in  length  (for  details  see 
Fig.  2-3).  Later,  8  more  hydrophones  were  deployed  for  a  total  of  32  in  the  horizontal 
array  (HA),  all  at  a  depth  of  60  m.  The  labels  near  some  hydrophones  on  Fig.  2-3 
and  2-2  indicate  the  nominal  direction  and  range  in  meters  of  the  sensors  as  surveyed 
relative  to  the  apex  of  the  site.  Two  USRD  J9  acoustic  sources^  were  deployed  and 
each  driven  with  a  dual  tone  signal  for  use  in  localizing  the  VLA  channels. 

A  subset  of  the  recorded  data  was  processed  in  real  time  to  continuously  update 
seismicity  displays.  The  data  were  passed  through  a  real-time,  focused  beamformer. 


^USRD  stands  for  “Underwater  Sound  Reference  Detachment”,  an  organization  that  is  part  of 
the  Naval  Research  Lab.  USRD  owns  various  pieces  of  acoustic  equipment  that  they  can  lend 
or  rent  to  other  institutions  for  temporary  use.  These  J9  sources  were  rented  from  them  by  the 
MIT/WHOI  group  for  the  SIMI-94  experiment. 
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Arrays  and  Sources 


Figure  2-2:  Seismo-acoustic  arrays  deployed  by  MIT/WHOI  on  and  around  the  SIMI- 
94  East  Camp  floe.  A  total  of  32  hydrophones  (shown  on  the  figure  by  the  circles 
with  crosses)  were  deployed  in  the  horizontal  array  in  a  combined  circular/cross  con¬ 
figuration.  The  labels  near  some  hydrophones  nominally  indicate  the  direction  and 
range  in  meters  of  the  sensors  as  surveyed  relative  to  the  array  apex  (the  origin  of 
the  coordinate  system  on  the  plot).  Only  the  names  of  the  outermost  hydrophones 
are  shown  on  this  figure  to  avoid  cluttering  it  up.  A  32  element  vertical  array  was 
deployed  at  the  horizontal  array  apex.  Two  J-9  acoustic  sources  (J9A  and  J9B)  were 
deployed  for  sensor  tracking.  Markers  with  labels  Ant  and  Gen  show  the  locations 
of  the  antenna  hut  and  generator  hut,  respectively.  The  location  of  the  geophone 
array  near  ice  mountain  called  Mount  Odyssey  is  shown  by  square  marker  with  label 
RLAM  near  it.  The  locations  of  the  ice  mountains  are  shown  by  triangle  markers. 
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Figure  2-3:  Horizontal  surveillance  hydrophone  array  of  the  SIMI-94  East  Camp.  The 
hydrophone  locations  are  shown  on  the  figure  by  the  circles  with  crosses.  The  labels 
near  some  hydrophones  indicate  the  nominal  direction  and  range  in  meters  of  the 
sensors  as  surveyed  relative  to  the  array  apex  (the  origin  of  the  coordinate  system  on 
the  plot).  A  32  element  vertical  array  was  deployed  at  the  horizontal  array  apex  (the 
label  VLA  reflects  this  fact). 


31 


Meters  North 


1.5e+03 


le+03- 


, . i . ^  ;•  ;Vi.*:  *i 


g'^'-  .i:!  ' . 

Jv\. ; . ^ . 

. ; r. .  JP^:.  ^  ?  ..  r.\  ^J.U..;,  .; 

. -i . V-;,-  •  _ 


0-^ 


>* 


•••. 


-le+034 


'  ••'V  ■»>»  • .  •  V.  t/t 

.  V»  i-:!  i  : 

':  ’i  :  .-v  . ;• . 0":- ;^.  v  r; . . :  - 

. . .•:-.;,..V,‘,  ,’*•  .;  ^g.^N  ’.'''^''T7'  ii; 

"•  :..‘[j,,,  •:  ..:  ,  V  -y  .,. 

,  ■  '*  ••  ■  ^'•:i  ^ . . . . . 

. . .  .:;:.4.i-.iSi.j 


-L5e+034 


-1.5e+03  -le+03 


T 

0 


|•>-14  db 
•>-20  db 
•>-26  db 
•>-32db 

•  >-38  db 
•>-44db 

•  >-50  db 


Meters  East 


le+03  1.5e+03 


Figure  2-4:  Ice  events  located  by  the  real-time  beamforming  on  horizontal  hydrophone 
array  between  098:23. OOZ  and  098:23. 59Z,  including  major  events  associated  with 
existing  crack  immediately  north  of  the  new  shear  ridge  (Mount  Odyssey  on  Fig.  2- 
2).  Such  real-time  beamforming  was  used  during  the  experiment  to  identify  the  active 
areas  on  ice. 


32 


developed  specifically  for  ice  event  detection  and  localization  by  Edward  K.  Scheer, 
and  hourly  event  maps  were  generated.  As  a  characteristic  example,  Fig.  2-4  shows 
events  detected  between  098:23.00Z  and  098:23.59Z,  including  major  events  associated 
with  an  existing  crack  immediately  N  of  a  new  shear  ridge  (Mount  Odyssey).  The 
size  of  the  markers  indicate  the  estimated  source  levels  of  the  individual  events.  As  a 
result  of  this  detection,  one  of  the  Radio  LAN  Acquisition  Module  (RLAM)  geophone 
clusters  was  deployed  in  the  vicinity  of  the  crack  and  recorded  the  continued  seismicity 
for  several  days  following.  This  example  is  introduced  only  to  show  how  active  areas 
were  identified  during  experiment  and  is  not  a  part  of  the  thesis  research. 

These  results  were  used  to  select  sites  for  deploying  autonomous  RLAMs  that 
acquired  geophone  data  in  the  near-field  of  apparently  active  sites.  This  data  is  the 
primary  focus  of  this  thesis. 

2.2.2  Geophone  Clusters 

Ice  motion,  either  as  a  result  of  natural  activity  or  induced  events,  was  measured  di¬ 
rectly  by  portable  RLAM  systems,  each  equipped  with  5  geophones  and  1  hydrophone. 
The  geophones  were  standard  Teledyne  units  mounted  in  an  OBS^  housing.  The  low 
frequency  cutoff  of  each  geophone  was  4.5  Hz  and  its  high  frequency  cutoff  was  above 
10  kHz^.  At  the  low  frequency  cutoff  the  geophones  were  damped  such  that  at  4.5  Hz 
the  amplitude  response  was  down  by  3  dB  and  continued  to  roll  off  as  frequency 
decreased  by  6  dB/octave.  Each  geophone  was  actually  a  three-component  cluster 
which  measured  three  orthogonal  components  of  particle  velocity  of  the  ice  motion. 
The  geophones  were  mounted  directly  on  the  hard  ice  using  slush,  which  immediately 
froze  the  units  to  the  ice.  In  case  of  snow  cover  present,  the  ice  was  first  cleared.  The 
typical  snow  cover  was  around  30  cm. 

These  systems  were  deployed  at  8  sites  in  the  vicinity  of  the  SIMI-94  East  Camp 
for  varying  times  throughout  the  experiment.  Data  was  telemetered  via  a  radio  local 

^  OBS  stands  for  Ocean  Bottom  Seismometer 

^  The  precise  value  of  the  high  frequency  cutoff  is  not  important,  because  the  maximum  sampling 
frequency  for  the  geophone  arrays,  as  far  as  I  know,  was  around  4  kHz. 
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Figure  2-5:  Directions  of  the  geophone  axes.  X  axis  is  oriented  along  the  main  axis 
of  the  geophone.  Both  X  and  Y  axes  lie  in  the  horizontal  plane.  Z-axis  is  the  vertical 
axis  with  positive  direction  upwards. 


area  network  to  a  receiving  system  located  at  the  MIT/WHOI  science  hut  at  the 
SIMI-94  East  camp  where  it  was  stored  on  8  mm  tape  in  a  format  similar  to  the 
‘surveillance’  system.  The  data  bandwidth  varied  from  400  to  1600  Hz  depending  on 
whether  1,  2  or  3  RLAM  units  were  operating  simultaneously. 

Data  was  acquired  during  selected  times  from  7  April  through  18  April  depend¬ 
ing  on  the  schedules  of  other  experiments,  such  as  induced  fracturing  and  ice  block 
drops  as  well  as  the  detection  of  natural  events.  During  the  remainder  of  the  experi¬ 
ment  through  23  April,  3  RLAM  units  were  kept  in  continuous  operation  near  shear 
ridge  ‘Mt.  Odyssey’,  the  remote  APL/UW^  site,  the  large  lead  to  the  north  of  the 
camp  (so  called  “ice  river”  site),  and  on  a  small  floe  involved  in  active  ridge  building 
approximately  4  km  east  of  the  camp  (so  called  “ice  island”  site). 


^  APL/UW  stands  for  Applied  Physics  Laboratory  of  University  of  Washington 
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2.3  Preparation  of  the  SIMI  data 


In  this  section  the  preparation  of  the  SIMI  geophone  data  is  discussed.  First,  the 
preprocessing  of  the  data  needed  to  obtain  absolute  levels  of  the  velocity  in  ice  is 
discussed.  After  that,  the  procedure  for  event  detection  is  presented.  Then  follows 
the  description  of  the  standard  procedure  for  event  localization.  The  method  begins 
with  the  estimation  of  time  delays  between  different  sensors,  followed  by  the  proce¬ 
dure  for  a  least-square  fit  to  the  interarrival  data  assuming  an  unknown  propagation 
speed  (the  section  in  text  labelled  “Search  in  Cartesian  coordinate  space”).  In  this 
procedure,  all  possible  locations  for  the  source  are  tested  and  the  choice  is  made  based 
on  the  least  deviation  between  calculated  and  experimental  interarrival  times.  Some¬ 
times,  it  is  advantageous  to  use  different  procedure  when  for  each  trial  propagation 
speed  (as  opposed  to  trial  source  location  in  previous  method)  the  standard  devia¬ 
tion  between  calculated  and  experimental  interarrival  times  is  determined  and  the 
value  of  propagation  speed  (and  corresponding  source  location)  which  minimizes  this 
difference  is  chosen.  The  description  of  that  method  (called  “Search  in  propagation 
speed  space”)  finishes  this  section. 

2.3.1  Preprocessing  of  geophone  data 

To  determine  the  actual  geophone  positions,  the  records  of  the  distances  and  an¬ 
gles  between  geophones  were  used.  These  records  were  processed  using  the  program 
SURVEY,  written  by  Prof.  H.  Schmidt  in  FORTRAN  programming  language.  This  pro¬ 
gram  used  a  simulated  annealing  method  to  determine  the  locations  of  the  geophones. 

Every  attempt  was  made  to  orient  the  main  axis  of  each  geophone  (axis  X  on 
Fig.  2-5)  along  the  direction  to  magnetic  North.  The  small  deviations  in  direction 
from  magnetic  North  at  each  geophone  and  magnetic  deviation  at  the  experimental 
site  were  accounted  for  using  conceptual  rotation  of  horizontal  components  of  velocity. 
So  if  the  original  (recorded)  horizontal  components  on  the  geophone  were  and 
y/gc  and  rotated  -  X*  (true  North  component)  and  F*  (true  West  component),  then 
the  conceptual  rotation  can  be  expressed  as 
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(2.1) 

yi 

sin(cVniag  “1“  “1“  ^ec^^^C^niag  “1” 

(2.2) 

where:  ai  -  deviation  of  main  axis  of  the  geophone  from  magnetic  North;  Ofmag  - 
magnetic  deviation  at  the  experimental  site. 

The  conversion  of  the  raw  data  to  true  units  (m/sec)  on  the  magnetic  tape 
involves  three  steps.  First,  the  RLAMIZE  routine  written  in  C-language  by  Edward  K. 
Scheer  converts  the  raw  data  value  to  volts  at  the  input  of  the  RLAM  analog-to-digital 
converters  but  before  fixed  gain  was  applied: 

void  rlamize (signed  short  *buf in, float  *bufout,int  npts) 

■C 

int  i; 

unsigned  short  exponent , one ; 

one  =  1; 

f or (i=0 ; i<npts ; i++) 

exponent  =  one«((bufin[i]&0x03)*3) ; 

buf out  [i]  =  (float)  (buf in  [i]  »2)  *4 . 5/  (8192 .  ♦  (float )  exponent* .  625)  ; 

} 

} 


This  routine  has  an  input  in  buf  in  and  output  in  buf  out  and  it  processes  npts  of 
2-byte  long  data  values.  To  correct  the  data  for  a  fixed  gain,  the  output  of  this  routine 
is  divided  by  factor  of  100.  The  last  step  in  obtaining  values  in  true  units  is  using  the 
geophone  sensitivity,  which  was  642  mV/in/sec,  or  25.3  V/m/sec  to  convert  volts  to 
m/sec. 
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2.3.2  Event  detection  procedure 

The  first  step  in  the  detection  algorithm  is  to  read  the  selected  block  of  data  B  (say 
20  sec.  long)  from  the  input  data  file  or  directly  from  the  magnetic  tape.  Then  the 
long-time  average  (LTA)  rms  energy  is  computed  for  each  geophone  and  for  each 
component  of  the  velocity.  As  a  result,  we  obtain  3  sets  of  N  values 
where  index  i  is  a  geophone  number  (i  =  1, . . .  ,  AT).  For  example, 


Lxi  — 


1  " 

i=i 


where: 

to  -  the  starting  time  of  block  B, 

M  -  the  number  of  samples  in  the  block, 
dt  -  the  sampling  interval  of  the  recording, 

Vxiih)  ~  the  value  of  x-component  of  the  velocity  of  the  ice  particle  recorded 
by  the  z-th  geophone  at  the  time  moment  ti, 

MLci  -  the  mean  of  the  x-component  of  velocity  recorded  by  z-th  geophone  over 
the  whole  block  B. 

Mri  is  defined  as 

1  ^ 

+  U  -  l)dt). 

j=i 

Single-channel  LTAs  are  needed  to  determine  the  Signal-to-Noise  Ratio  (further  re¬ 
ferred  to  as  SNR).  Note  that  this  allows  adaptation  to  varying  noise  levels  in  time  as 
we  move  through  the  data  set  by  reading  a  new  block  from  the  file  or  directly  from 
the  tape. 

Averaging  LTAs  also  across  the  geophones,  we  obtain  now  3  values:  L^,  Lj,, 
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.  So,  for  example, 


1  ^ 

i=l 


where  N  is  the  number  of  the  geophones. 

In  addition,  we  keep  track  of  long-time  means  averaged  across  the  geophones 
Mr ,  IVl^ ,  Mr .  For  example, 


1  TV  M 

i=l  j=l 

where  to  “  the  starting  time  of  block  B,  M  -  the  number  of  samples  in  the  block,  dt  - 
the  sampling  interval  of  the  recording,  and  Va,i{ti)  -  the  value  of  ar-component  of  the 
velocity  of  the  ice  particle  recorded  by  the  t-th  geophone  at  the  time  moment  ti . 

We  choose  a  window  length  of  analysis  W  such  that  it  is  long  enough  to 
encompass  the  typical  signal  lengths  expected  in  the  recorded  time  series.  The  other 
requirement  for  the  length  of  window  is  that  it  should  allow  a  sufficient  time  for  the 
signal  to  appear  on  all  the  channels  in  the  analysis  window. 

Next,  the  short-time  average  (STA)  rms  energy  is  computed  for  each  geophone 
and  each  component  of  the  velocity.  As  a  result  we  obtain  3  sets  of  N  values  S^i,  Syi, 
Sxi,  where  index  f  is  a  geophone  number  (i  =  1, . . .  ,  N). 

Using  per-channel  STAs  and  averaged  over  the  channels  LTAs  (Lj,  L^,  L^) 
we  determine  whether  the  particular  channel  (particular  component  of  the  particular 
geophone)  is  “noisy”  or  not.  For  example,  let  us  consider  the  ar-component  of  velocity 
recorded  by  the  i-th  geophone.  Then,  if 

(Mr  —  21Lj;)  <  Sxi  <  (Ma;  -b  21*^)  , 

this  particular  channel  is  declared  “non-noisy”  and  can  be  used  for  detection  purposes 
in  window  W.  Otherwise,  it  is  discarded  for  detection  in  this  time  window. 

Next,  we  determine  whether  the  data  within  the  analysis  window  is  signal 
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or  noise  only  for  “non-noisy”  channels  (see  above).  This  is  done  by  computing  the 
SNR  for  each  receiver  within  the  analysis  window.  Basically,  this  SNR  is  nothing  more 
than  the  ratio  of  the  short-time  average  (STA)  to  long-time  average  over  the  particular 
recorded  channel  of  the  data.  For  example,  for  the  x-component  of  z-th  geophone, 
this  ratio  will  be  SxilL^i,  where  S^i  is  STA  over  window  W  for  x-component  of  z-th 
geophone,  and  Lxi  is  LTA  over  block  B  for  x-component  of  z-th  geophone.  Only  if 
such  defined  SNR  is  greater  than  a  preset  threshold,  the  channel  registers  that  the 
signal  is  present  in  the  current  time  window. 

A  variation  on  this  technique  uses  as  the  noise  level  the  STA  over  the  window 
with  minimal  energy.  Such  a  window  is  found  by  scanning  the  whole  block  for  the 
“quietest”  window  TTmin  without  regard  to  the  particular  velocity  component  (that 
means  one  window  for  all  3N  channels).  In  that  case,  the  SNR  in  window  W  for 
the  above  example  defined  as  Sxi/Smin,  where  ^min  is  the  STA  over  the  window  with 
minimal  energy  Wmm- 

A  detection  is  declared  when  at  least  half  of  the  geophones  register  a  signal 
in  the  window  W  for  at  least  one  of  the  components  of  the  velocity  (X-,  Y—  or  Z- 
component).  If  that  is  not  the  case,  the  recorded  time  series  is  considered  to  be  noise, 
and  therefore  we  move  to  the  next  analysis  window. 

2.3.3  Event  localization  scheme 

Time  delay  estimation 

For  each  event,  time  delays  of  the  event  arrival  at  different  geophones  are  calculated 
using  a  standard  cross-correlation  procedure.  A  reference  channel  is  first  selected, 
usually  the  one  where  the  event  arrives  first.  The  signal  in  that  channel  is  then  cross- 
correlated  with  the  corresponding  signals  on  all  the  other  geophones  on  which  the 
event  has  been  detected.  The  time  delay  is  the  time  at  which  the  cross-correlation 
coefficient  attains  an  absolute  maximum. 
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Figure  2-6:  Simple  receiver  configuration  for  establishing  the  relations  between  arrival 
times.  Geophones  are  located  at  the  points  A  and  B  and  the  source  at  point  S.  The 
distance  between  geophones  is  Rab-  The  distance  from  geophone  A  to  the  source  is 
Ra,  while  the  distance  from  geophone  B  to  the  source  is  Rb- 


The  cross-correlation  function  is  defined  as 


Kxi ,xi {r)  =  E [xi{t)xj {t  -  r)] ,  (2.3) 

where  E  is  an  expectation  operator.  Due  to  finite  duration  of  the  time  series,  only 
an  estimate  of  K  can  be  obtained.  The  value  of  r  that  maximizes  the  expression  in 
Equation  (2.3)  is  the  time-delay  (or  arrival  time)  estimate. 

There  are  some  intrinsic  checks  on  the  value  of  delays  between  the  different 
sensors.  Let  us  consider  two  sensors  (geophones)  located  at  points  A  and  B  and  the 
source  at  point  S  (see  Fig.  2-6).  The  propagation  speed  is  assumed  to  be  cq.  Then, 
according  to  well-known  geometrical  inequalities,  the  distance  between  the  geophones 
Rab  is  not  greater  than  the  sum  of  distances  from  each  geophone  to  the  source  (further 
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Figure  2-7:  Two  examples  of  the  cross-correlation  between  X-components  of  the  geo¬ 
phones.  At  the  top  plot  the  original  time  series  start  at  37.5  sec.  from  the  beginning 
of  the  tape  RLAM-26,  on  the  bottom  one  -  at  48  sec.  from  the  beginning  of  the  tape 
RLAM-26.  In  both  cases  the  cross-correlations  were  carried  out  with  geophone  #2 
being  a  reference  channel. 
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referred  to  as  Ra  and  Rb)'- 


Rab  <  Ra  +  Rb- 

On  the  other  hand,  the  absolute  value  of  the  difference  between  the  last  two  ranges 
is  not  greater  than  Rab'- 


\Ra  —  Rb\  <  -Rab- 


Dividing  both  sides  of  the  last  inequality  by  Cq,  we  get: 

\U  -  (bI  <  — , 

where  Ia  and  are  propagation  times  from  the  source  to  the  corresponding  receiver. 
The  distances  between  the  geophones  are  on  the  order  of  50  m.  For  the  typical  values 
of  compressional  and  shear  wave  propagation  speeds  in  ice  (3500  m/s  and  1800  m/s), 
we  obtain  the  following  maximum  allowed  delays  between  geophones: 

Ate  ~  14  msec.  (2.4) 

Atg  ~  28  msec.,  (2.5) 

where  Ate  is  the  maximum  allowed  delay  between  geophones  for  the  compressional 
waves,  and  Atg  is  the  maximum  allowed  delay  between  geophones  for  the  shear  waves. 

Two  examples  of  the  cross-correlation  are  shown  on  Fig.  2-7.  In  both  examples 
the  outputs  of  X-axes  of  the  geophones  were  used  and  geophone  #2  was  chosen  as  a 
reference  channel. 

Least-square  fit  for  interarrival  times:  Search  in  Cartesian  coordinate 
space 

In  this  method  a  rectangular  grid  is  constructed,  with  cell  size  10  m  x  10  m.  The 
source  is  placed  at  each  grid  point.  The  distance  from  the  source  location  to  each 
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Figure  2-8:  Example  of  the  least-square  fit  to  arrival  times  on  4  geophones.  The 
measured  arrival  times  are  plotted  versus  ranges  to  the  best-fit  source  location.  The 
location  of  event  was  estimated  to  be  110  m  to  the  east  and  430  m  to  the  north  from 
the  array  apex.  Standard  deviation  -  0.0022.  Estimated  wave  speed  -  3300  m/s. 
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Figure  2-9:  Locations  of  events  determined  by  standard  least-square  fit.  Triangles 
mark  geophone  locations,  filled  circles  —  event  locations.  The  numbers  nearby  markers 
show  the  moment  of  time  (from  the  beginning  of  the  RLAM  tape  26)  at  which  event 
occured  (rounded  to  seconds). 


geophone  is  computed,  as  well  as  corresponding  time  delays.  The  standard  deviation 
from  the  measured  delay  values  is  then  calculated  as 


(2.6) 


where  N  is  the  number  of  the  geophones  and  Ti  and  are  the  measured  and 
estimated  time  delays,  respectively.  The  propagation  speed  is  also  estimated  as  the 
slope  of  the  best-fit  regression  line  through  the  points  with  Ri  being  the 

range  for  each  geophone  from  the  source  location.  The  parameters  A  and  B  of  the 
least-square  fit  Ti  =  ARi  -f  B,  can  be  determined  as 


B  = 


N  N  N 

J2nRi  - 

t=l _ i=l  t=l 

N  /  N  \‘^ 

JlRl-iiERi) 

i=l  \i=l  / 

N  N 

fin  -  AflRi 

i=l _ i=l 

N 


(2.7) 


(2.8) 


Then  for  the  propagation  speed  cq,  we  have  cq  =  1/A,  and  the  Equation  (2.6)  can  be 
rewritten  as 


CX'T  — 


E."  1  (n  -ARi-  Bf 
N 


For  the  example  of  least-square  fit  shown  on  Fig.  2-8,  the  value  for  the  wave 
propagation  speed  was  estimated  to  be  3300  m/s,  which  is  realistic  for  typical  Arctic 
ice  conditions. 

Several  events  located  by  this  method  at  the  beginning  of  the  REAM  tape 
26  are  shown  on  the  Fig.  2-9.  We  see  from  the  figure  that  events  occurring  at  the 
moments  of  time  168  sec,  186  sec  and  214  sec  after  the  beginning  of  the  tape  appear 
to  lie  on  the  same  arc,  almost  at  the  limit  of  the  region  of  reliable  localization  by  the 
geophone  array. 
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Other  events  at  the  beginning  of  the  same  tape  I  grouped  separately  because 
their  waveforms  had  similar  features  which  repeated  from  one  event  to  another  (see 
Fig.  2-10).  Compare,  for  example  the  degradation  of  the  shape  of  the  arrival  on  the 
first  geophone  as  compared  to  other  geophones  for  all  shown  events.  The  locations 
of  events  corresponding  to  these  time  series  (including  one  not  shown  on  Fig.  2-10) 
are  shown  on  Fig.  2-11.  I  call  this  group  of  5  events  Cluster  #1.  In  this  group  all 
events  are  at  approximately  the  same  azimuth  and  their  locations  are  quite  close  to 
each  other  (except  the  last  one). 

Least-square  fit  for  interarrival  times:  Search  in  propagation  speed  space 

The  main  difference  from  the  previous  method  is  that  the  search  is  carried  out  in 
propagation  speed  space  as  opposed  to  the  search  in  Cartesian  coordinate  space  for 
trial  source  locations.  In  this  method  we  try  different  values  of  propagation  speed 
and  again  minimize  the  standard  deviation  over  the  search  space.  For  each  trial  value 
of  propagation  speed  cq  in  the  ice  plate,  the  arrival  time  ti  at  the  z-th  geophone  can 
be  expressed  as  [20] 

ti  =  ts  +  —\/{Xs-  XiY  +  {ys  -  yiY  -1-  {zs  -  ZiY,  (2.9) 

where  ts  is  the  actual  moment  of  time  when  the  ice  event  located  at  {xg,  ys,  Zg)  hap¬ 
pened.  The  N  geophones  are  assumed  to  be  positioned  at  (xi.yi.Zi)  with  (i  = 
1, ...  ,N).  Usually,  the  geophones  are  deployed  approximately  at  the  same  depth 
(further  referred  to  as  2:0  )• 

Multiplying  both  sides  of  Equation  (2.9)  by  cq,  squaring  them,  and  then  sub¬ 
tracting  from  the  equations  for  (i  =  2, . . .  ,  N)  the  equation  for  first  geophone  (i  =  1), 
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Figure  2-10:  Stacked  time  series  of  the  X-component  (North)  of  the  velocity  for  four 
events  at  the  beginning  of  the  RLAM  tape  26.  The  time  series  corresponding  to 
the  geophone  #1  are  at  the  top.  All  events  have  strikingly  similar  features  in  their 
waveforms.  Compare  the  degradation  of  the  shape  of  the  arrival  on  the  first  geophone 
as  compared  to  others  for  all  shown  events. 
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Figure  2-11:  Locations  of  events  time  series  for  which  have  striking  similarities  (com¬ 
pare  Fig.  2-10).  Triangles  mark  geophone  locations,  filled  circles  -  event  locations. 
The  numbers  nearby  markers  show  the  moment  of  time  (from  the  beginning  of  the 
RLAM  tape  26)  at  which  event  occured  (rounded  to  seconds).  In  this  group  all  events 
are  at  approximately  the  same  azimuth  and  their  locations  are  quite  close  to  each 
other  (except  the  last  one). 
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we  obtain  the  following  system  of  {N  —  1)  equations: 


2(x2  -  a;i)  2(y2  -  Vi)  -2<^(i2  -  <i) 

f  > 

Xg 

{xl  -\-yl-  (^tl)  -  sf 

:  :  : 

< 

Vs 

^  =  < 

: 

2{xif  -  xi)  2(yjv  -  yi)  -2c^(%  —  ti) 

where 


Si 


=  \[A 


+  yi 


c^tl 


So  we  need  to  solve  this  system  of  {N  —  1)  equations  for  three  unknowns  Xg,  ys  and 
tg.  In  matrix  notation  the  above  system  takes  the  form  of 

AX  =  B,  (2.10) 

where  matrix  A  contains  all  coefficients  from  the  lefthand  side  of  the  system  of  equa- 
tions,  the  vector  X  contains  all  unknowns  and  the  vector  B  contains  all  coefficients 
from  the  righthand  side  of  the  original  system  above. 

The  solution  of  Equation  (2.10)  can  be  obtained  using  the  least-square  solution 
described  in  Appendix  A: 


X  =  (A’’A)"‘  (a’’b)  .  (2.11) 

Note  that  this  type  of  solution  works  only  when  the  matrix  (A^A)  is  non-singular 
(i.e.  it  has  an  inverse).  Under  that  condition,  we  can  determine  the  corresponding 
source  location  {xg,yg)  and  the  standard  deviation  between  calculated  and  experi¬ 
mental  interarrival  times  (defined  similarly  to  Equation  (2.6))  for  each  value  of  trial 
propagation  speed  cq.  The  location  of  the  source  at  the  output  of  this  method  will 
be  the  location  corresponding  to  the  value  of  the  trial  propagation  speed  cq  which 
minimizes  the  standard  deviation  between  calculated  and  experimental  interarrival 
times. 
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2.4  Flexural  wave  arrival 


In  this  section  some  results  for  the  flexural  wave  arrival  will  be  discussed.  First,  the 
method  for  determining  dispersion  characteristics  of  the  arrival,  the  multiple  filter 
technique  (further  referred  to  as  MFlt),  will  be  described.  After  that,  the  disper¬ 
sion  pattern  of  the  flexural  wave  arrival  determined  by  MFlt  method  will  be  dis¬ 
cussed.  Finally,  the  sense  of  motion  in  flexural  waves  will  be  presented  based  on 
velocity /displacement  hodographs  of  motion  in  the  flexural  wave  arrival. 

2.4.1  Multiple  filter  technique 

The  first  step  in  the  multiple  filter  technique  by  Dziewonski  et  al.  [9]  is  to  pass  the 
data  through  a  bank  of  Gaussian  narrowband  filters  covering  the  frequency  range 
of  interest.  The  authors  use  Gaussian  filters  in  order  to  obtain  an  optimal  balance 
between  frequency  and  time  resolution.  The  general  expression  for  one  of  such  filters 
in  the  frequency  domain  is 


Hn{u))  =  e““(“)',  (2.12) 

where  is  the  center  frequency  of  the  n-th  filter  and  a  is  the  parameter  that  controls 
the  resolution  of  the  filter.  In  original  paper,  a  filter  bank  with  constant  relative 
bandwidth  was  used,  but  here  it  was  more  convenient  to  use  filters  with  constant 
bandwidth. 

In  the  next  stage  of  this  method,  the  analytic  signal  for  each  frequency  band 
used  is  calculated.  The  analytic  signal  is  defined  as: 

S{t]  iVc)  =s(t;  ojc)  +  /H{s(i;  ofc)}  (2.13) 

S{t-,co,)=A{t-uj,)e^‘l>^^’‘^-^\  (2.14) 


where: 

iOc  -  center  frequency  of  the  narrowband  filter. 
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s{t]  cOc)  -  filtered  signal, 

H  -  Hilbert  transform, 

S{t;u}c)  -  analytic  signal, 

A{t-,  Uc)  -  instantaneous  amplitude  of  the  signal, 
cJc)  -  instantaneous  phase  of  the  signal. 

Then,  the  group  arrival  time  for  the  center  frequency  of  each  band  is  estimated 
as  the  time  at  which  the  maximum  of  the  complex  envelope  (or  the  absolute  value  of 
the  analytic  signal)  occurs.  The  group  arrival  time  can  be  expressed  as 


where: 

R  -  distance  from  source  to  receiver, 

Vg  -  group  velocity. 

Using  the  above  expression,  we  can  now  estimate  the  group  velocity  of  the 

signal. 

2.4.2  Flexural  wave  dispersion 

Fig.  2-12  shows  the  location  of  the  event  that  produced  vertical  velocities  charac¬ 
teristic  of  a  flexural  wave  dispersion  pattern  (first  arrive  higher  frequencies,  later  - 
lower  frequencies).  The  time  series  of  vertical  velocities  for  this  event  are  shown  on 
Fig.  2-13.  These  times  series  were  processed  using  the  MFlt  technique.  The  result  of 
this  processing  for  the  geophone  #5  is  shown  on  Fig.  2-14.  On  this  plot  the  arrival 
times  were  converted  to  group  velocities  using  the  distance  from  the  geophone  to  the 
event  location  and  the  time  the  event  occured,  obtained  previously. 

The  general  trend  can  be  seen  even  on  this  individual  plot,  but  quantitative 
results  are  better  obtained  from  averaging  across  geophones.  This  is  done  in  the 
following  way:  the  arrival  time  corresponding  to  the  maximum  spectral  amplitude  at 
each  frequency  for  each  geophone  separately  is  determined,  this  time  is  converted  to 
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Figure  2-12:  Location  of  the  event  that  produced  the  vertical  velocities  characteristic 
for  the  flexural  wave  dispersion  pattern  (flrst  arrive  higher  frequencies,  later  -  lower 
frequencies).  Event  location  is  shown  by  filled  circle  marker,  while  geophone  loca¬ 
tions  are  shown  by  unfilled  triangles.  This  event  location  was  determined  using  the 
least-square  fit  for  the  arrival  times  determined  by  cross-correlation  of  the  vertical 
component  of  the  velocity  registered  by  the  geophones. 
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Figure  2-13;  Stacked  time  series  of  vertical  velocities  for  the  event  that  produced  the 
vertical  velocities  characteristic  for  the  flexural  wave  dispersion  pattern  (first  arrive 
higher  frequencies,  later  -  lower  frequencies).  The  time  series  recorded  by  the  first 
geophone  are  at  the  bottom,  while  time  series  corresponding  to  the  last  geophone  - 
at  the  top. 
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Figure  2-15:  Dispersion  curves  of  flexural  wave.  Five  different  type  of  markers  cor¬ 
respond  to  different  geophones  (see  legend).  Average  across  geophones  is  shown  on 
this  plot  as  a  solid  line.  Note,  that  these  results  should  be  believed  only  till  about 
20Hz,  because  higher  frequencies  are  absent  in  time  series,  as  it  can  be  seen  from  the 
Fig.  2-13. 

the  group  velocity  using  the  known  range  to  source,  then  5  different  curves  (different 
markers  corresponding  to  different  geophones)  are  plotted  as  group  velocity  versus 
frequency  (see  Fig.  2-15).  The  average  across  geophones  is  shown  on  this  plot  as  a 
solid  line  (in  one  case,  at  18  Hz,  I  removed  one  outlier).  This  dispersion  curve  should 
be  believed  only  until  about  20  Hz,  because  higher  frequencies  are  absent  in  time 
series,  as  can  be  seen  from  the  Fig.  2-13. 
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Figure  2-16;  Prograde  and  retrograde  senses  of  the  particle  motion.  According  to 
[1],  the  convention  for  distinguishing  between  prograde  and  retrograde  motion  can 
perhaps  best  be  remembered  in  terms  of  a  ball  rolling  along  the  top  of  the  half-space 
from  source  to  receiver.  The  sense  of  rotation  of  the  ball  is  prograde,  while  for  the 
Rayleigh  wave  at  the  surface  the  particle  motion  is  retrograde. 

2.4.3  Sense  of  motion  in  flexural  wave 

The  different  modes  of  particle  motion  can  be  distinguished  by  their  so  called  “sense  of 
particle  motion” .  According  to  [1],  the  convention  for  distinguishing  between  prograde 
and  retrograde  particle  motion  can  perhaps  best  be  remembered  in  terms  of  a  ball 
rolling  along  the  top  of  the  half-space  from  source  to  receiver  (see  Fig.  2-16).  When  the 
wave  is  propagating  to  the  right,  prograde  sense  of  motion  corresponds  to  clockwise 
motion,  while  retrograde  sense  -  counterclockwise.  For  the  fundamental  mode  of  the 
Rayleigh  wave  at  the  surface,  the  sense  of  particle  motion  is  retrograde. 

Fig.  2-17  shows  the  hodographs  of  velocity  and  displacement  on  one  of  the 
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Figure  2-17:  Hodographs  of  velocity  and  displacement  in  the  flexural  wave  arrival 
for  geophone  #5.  Top  plot  shows  hodograph  of  velocity,  while  the  bottom  plot  - 
the  hodograph  of  displacement.  Each  hodograph  shows  vertical  component  of  mo¬ 
tion  (velocity  or  displacement)  versus  radial  component  of  motion  (in  the  direction 
away  from  the  source).  The  circle  markers  correspond  to  the  start  of  the  time  in¬ 
terval  visualized  on  the  hodographs.  On  the  hodograph  of  velocity  (top  plot)  two 
additional  moments  in  time  are  marked  on  plot  by  filled  squares  to  resolve  ambiguity 
of  hodograph  behavior  near  its  first  self-intersection.  Numbers  near  markers  show 
approximate  times  in  msec  from  the  start  of  the  plotting  interval. 
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Figure  2-18:  Phase  shift  between  vertical  and  radial  components  of  motion  recorded 
by  the  geophone  #5.  Each  curve  shows  the  time  dependence  of  the  phase  shift  for 
different  frequency  (from  12  to  18  Hz).  Note,  that  all  the  curves  show  the  phase  shift 
close  to  —90°  at  time  interval  between  35.5  sec  and  35.7  sec.  This  phase  difference 
corresponds  to  the  retrograde  sense  of  motion,  because  vertical  component  is  behind 
of  radial  component  for  quarter  of  period  or,  in  other  words,  -90°. 

geophones  during  the  flexural  wave  arrival  discussed  above.  Both  hodographs  are  for 
motion  in  vertical  radial  plane.^  The  circle  markers  correspond  to  the  start  of  the 
time  interval  visualized  on  the  hodographs.  The  direction  of  the  wave  propagation 
is  shown  on  the  plot  with  horizontal  arrows  pointing  to  the  right.  The  arrows  on 
the  hodographs  show  the  sense  of  motion.  In  this  case,  it  appears  to  be  retrograde 
motion  as  for  the  fundamental  mode  of  a  flexural  wave  (the  only  mode  that  could 
possibly  exist  at  such  low  frequencies). 

The  phase  differences  between  vertical  and  radial  components  of  motion  on 


®  I  use  this  name  to  denote  the  vertical  plane  containing  both  the  receiver  and  the  source. 
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several  frequencies  confirm  the  retrograde  sense  of  motion  in  this  flexural  wave  arrival 
(see  Fig  2-18).  On  this  figure  the  different  lines  correspond  to  different  frequencies 
(from  12  to  18  Hz,  see  the  legend).  All  curves  on  this  plot  are  in  the  vicinity  of  —90° 
phase  shift  at  the  time  interval  between  35.5  sec  and  35.7  sec.  This  phase  difference 
corresponds  to  the  retrograde  sense  of  motion,  because  the  vertical  component  lags 
the  radial  component  by  90  degrees. 


59 


Chapter  3 


Polarization  processing 

3.1  Overview 

Polarization  methods  are  widely  used  in  studies  of  vector  wave  fields  to  exploit  the 
specific  relations  between  characteristics  of  wave  motion  in  orthogonal  directions.  In 
earthquake  seismology  and  geophysical  prospecting  polarization  processing  is  used 
to  separate  physically  different  wave  types  in  the  Earth  crust.  The  same  purpose 
could  be  achieved  for  ice  events  by  adapting  the  methods  of  polarization  processing 
from  earthquake  mechanics.  The  specific  method  chosen  for  adaptation,  so  called 
“Motion  Product  Detector”  (MPD)  method,  is  described  in  the  first  section  of  this 
chapter.  Two  schemes  of  data  processing  in  this  method  are  designed  to  isolate  waves 
with  predominantly  rectilinear  and  elliptical  type  of  particle  motion  respectively.  The 
details  of  the  implementation  of  the  MPD  method  for  ice  event  data  processing  are 
discussed  in  the  second  section.  As  one  of  the  byproducts  of  data  processing  in  this 
method  the  bearing  of  the  source  is  also  obtained.  Using  this  information  the  source 
location  can  also  be  determined  through  a  triangulation  procedure.  The  results  of 
the  discrimination  between  different  wave  types  and  event  localization  are  discussed 
in  the  next  section.  The  last  section  is  devoted  to  the  summary  and  conclusions  for 
the  polarization  processing  results. 
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Figure  3-1:  Stacked  time  series  of  the  horizontal  (X)  component  of  velocity  on  5 
geophones  which  serves  as  an  example  of  the  case  when  cross-correlation  of  time 
series  performs  poorly. 

3.2  Motion  Product  Detector  (MPD)  seismogram 
method 

As  it  was  already  mentioned,  one  of  the  first  stages  in  the  processing  of  the  record¬ 
ings  of  the  natural  ice  events  is  localization  and  identification  (or  discrimination)  of 
the  events.  The  traditional  way  to  do  localization  by  least-square  fit  to  the  arrival 
times  on  the  different  geophones  occasionally  fails,  probably,  due  to  overlapping  of 
the  different  wave  types  or  waves  generated  by  different  ice  events  (see,  for  example. 


61 


Polarization  plane 


S -  source 
R  -  receiver 


Figure  3-2:  Plane  of  polarization  of  signals  considered  in  MPD  processing.  This 
vertical  plane  contains  both  the  source  and  the  receiver. 


Fig.  3-1).  In  these  cases  it  is  sometimes  possible  to  achieve  localization  using  the  Mo¬ 
tion  Product  Detector  (MPD)  seismogram  method  inspired  by  J.E.  White’s  original 
work  [53].  This  method  allows  to  isolate  vertically  polarized  waves  (for  definition  of 
their  polarization  plane  see  Fig.  3-2)  with  two  kinds  of  particle  motion  trajectories: 
rectilinear  and  elliptical  (see  Fig.  3-3).  In  the  seismology  literature  the  term  “po¬ 
larization”  is  used  to  distinguish  between  these  two  kinds  of  particle  motion.  The 
difference  between  rectilinear  polarization  and  elliptic  polarization  is  the  phase  shift 
oi'K/2  between  the  vertical  and  horizontal  components  of  motion  in  the  latter  case, 
as  compared  to  no  phase  shift  at  all  for  rectilinear  polarization.  Usually,  body  waves 
have  rectilinear  polarization,  while  surface  waves  (e.g.  Rayleigh)  and  flexural  waves 
have  elliptic  polarization  [17],  but  due  to  interaction  with  the  free  surface  of  the  ice 
sheet  sometimes  SV-waves  acquire  elliptic  polarization.  So  the  ability  to  separate 
these  two  polarizations  would  allow  to  discriminate  between  different  wave  types. 


Rectyinear  polarization  Elliptical  polarization 


Figure  3-3:  Two  different  kinds  of  polarization  of  seismo-acoustic  signals  in  the  verti¬ 
cal  plane:  rectilinear  and  elliptical  motion.  The  first  one  is  characterized  by  in-phase 
horizontal  and  vertical  components  of  motion,  the  second  by  the  phase  shift  of  7r/2 
between  these  components. 

3.2.1  Procedure  for  isolating  rectilinear  polarization  in  data 
recordings  by  the  MPD  method 

To  separate  the  compressional  wave  each  of  the  horizontal  components  of  ice  motion 
registered  by  geophone  is  multiplied  by  the  vertical  component  and  averaged  by  inte¬ 
gration  over  a  certain  time  window  (see  Fig.  3-4a).  The  optimal  window  length  was 
determined  empirically  to  be  20  msec.  The  result  of  averaging  the  product  of  the 
horizontal  X-component  with  the  vertical  Z-component  is  called  Xhv?  while  the  re¬ 
sult  for  the  horizontal  F-component  is  called  Ihv-  To  visualize  the  processing  results 
Fhv  is  plotted  versus  Xhv-  As  a  result  of  the  time  averaging  we  obtain  an  enhanced 
output  from  the  compressional  wave  due  to  the  in-phase  horizontal  and  vertical  com¬ 
ponents  of  motion  produced  by  that  wave,  while  the  motion  produced  by  other  waves 
will  be  substantially  attenuated  due  to  the  different  phase  relationship  between  the 
vertical  and  horizontal  components  in  those  waves. 

So  in  the  case  when  the  compressional  wave  is  present,  the  output  of  this  HV 
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(a) 


Figure  3-4:  Schemes  of  processing  in  HV/HiV  MPD  methods:  (a)  HV  MPD 
method;  (b)  HiV  MPD  method:  the  block  with  %  in  it  corresponds  to  the  operation 
of  phase  shifting  vertical  component  for  7r/2  radians  which  is  achieved  using  Hilbert 
transform. 

processor  (name  suggested  by  J.E.  White  in  [53])  will  approximate  the  straight  line 
passing  through  the  origin  and  randomly  scattered  points  otherwise.  The  straight 
line  in  case  of  successful  detection  will  determine  the  bearing  of  the  source,  because 
we  consider  here  only  vertically  polarized  waves.  Recall  that  the  polarization  plane 
is  the  vertical  plane  that  contains  both  the  receiver  and  the  source  of  the  ice  event 
(see  Fig.  3-2).  As  a  result,  this  line  gives  also  the  projection  of  bearing  of  the  source 
on  the  horizontal  plane.  Using  the  output  of  such  a  processor  for  several  geophones 
we  can  easily  triangulate  the  source  location  in  the  horizontal  plane. 
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3.2.2  Procedure  for  isolating  elliptic  polarization  in  data  record¬ 
ings  by  the  MPD  method 

A  similar  idea  is  employed  in  the  HiV  processor^  (see  Fig.  3-4b).  The  main  difference 
between  HV  and  HiV  processors  is  that  in  the  latter  the  vertical  component  is 
shifted  in  phase  by  7r/2  radians.  Originally,  I  used  numerical  differentiation  to  achieve 
this  phase  shift,  but  due  to  the  nonuniform  frequency  response  of  the  numerical 
differentiation  operation,  I  later  switched  to  the  Hilbert  transform  to  achieve  that 
purpose.  The  result  of  averaging  the  product  of  the  horizontal  X-component  with 
the  inverted  vertical  Z-component  is  called  Xhivj  while  the  result  for  the  horizontal 
y-component  is  called  yniv-  We  can  express  them  at  the  time  t  as: 

t+T,/2 


r{t)=  j 


t-T,l2 

t+T,/2 

Yuvit)  =  J  VyV^dt 

t-Ti/2 

t+Ti/2 

Xmv{t)  =  J  Wx<^^dt 

t-T,/2 

t+T,/2 

Fkiv(^)  =  J  VyV^^^dt, 

t-Ti/2 

where  Vx,  Vy,  are  the  corresponding  components  of  the  velocity,  is  the  phase- 
shifted  vertical  component,  T/  is  the  integration  (or  averaging)  interval.  For  discrete 


HiV  according  to  [53]  stands  for  Horizontal  and  inverted  Vertical  components. 
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signals  the  above  relations  transform  to 


M/2 

Vx{tj  +  iAt)vz{tj  +  iAt)  (3.1) 

i=-M/2 

M/2 

Yjiv{tj)=  ^  Vy{tj  +  iAt)v^{tj  +  iAt)  (3.2) 

i=-M/2 

M/2 

-^Hiv(0  =  ^2,  (3-3) 

i--M/2 

M/2 

^Hiv(0  =  ^2  Y  iAt)v^^^{tj  +  lAt)  (3.4) 

i=-M/2 

To  visualize  the  processing  results,  Tniv  is  usually  plotted  versus  Xhiv- 

As  a  result  at  the  output  of  such  a  processor  we  would  get  a  greatly  enhanced 
SV-wave  (or  Rayleigh  wave),  because  the  phase-shifted  vertical  component  will  now 
be  in  phase  with  horizontal  components.  The  other  waves  will  produce  a  substan¬ 
tially  attenuated  output  due  to  the  different  phase  relationship  between  vertical  and 
horizontal  components  in  those  waves. 

So  as  a  result  the  HV  MPD  method  can  be  used  as  compressional  wave  de¬ 
tector,  while  its  HiV  counterpart  is  used  for  SV-waves  or  Rayleigh  waves.  Especially 
useful  is  the  ability  of  this  method  to  achieve  simultaneously  both  the  separation  of 
the  different  polarizations  of  seismic  waves  and  the  bearing  of  the  source  of  the  ice 
event. 


3.3  Data  processing  scheme 

Using  the  MPD  polarization  processing  method  the  development  of  ice  fractures  in 
the  array  near-field  was  successfully  tracked  in  the  time  and  spatial  domains. 

The  visualization  procedure  is  a  little  bit  different  from  the  one  described  for 
both  processors.  The  reason  for  this  change  is  the  presence  of  two  coordinate  systems 
in  this  particular  experimental  setup:  the  first  one  (further  referred  to  as  So)  is  the 
Cartesian  coordinate  system  in  which  the  locations  of  sensors  (geophones)  are  defined, 
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Array  coordinate  system 


-40  -20  0  20  40 

Meters  East 


G^ophone  coordinate  system 


Figure  3-5:  Two  coordinate  systems  relevant  to  the  geophone  data.  The  top  plot 
shows  the  coordinate  system  of  the  geophone  array  (referred  to  as  So  in  the  text). 
Markers  show  the  locations  of  geophones.  The  origin  is  usually  the  location  of  the 
RLAM  antenna.  The  direction  of  the  x-axis  is  due  East,  that  of  the  y-axis  is  due 
North.  The  bottom  plot  shows  the  plan  view  of  a  geophone  with  its  internal  coordinate 
system  (this  system  is  called  Sg  in  the  text).  The  direction  of  the  x-axis  (main 
geophone  axis)  is  due  North,  that  of  the  y-axis  is  due  West.  These  directions  are 
identical  for  all  geophones  due  to  the  “conceptual”  rotation  procedure. 


67 


Figure  3-6:  Processing  results  of  the  Event  F5-3  (see  in  the  text)  with  the  HV  MPD 
method.  For  definitions  of  Xmpd  and  Fmpd  see  Equations  (3.5)  and  (3.6).  Different 
markers  correspond  to  different  geophones  (see  legend).  For  correspondence  between 
geophone  names  and  geophone  numbers  used  throughout  the  text  see  Fig.  3-8. 
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Figure  3-7:  Processing  results  of  the  Event  F5-3  (see  in  the  text)  with  HiV  MPD 
method.  For  definitions  of  Xmpd  and  Impd  see  Equations  (3.7)  and  (3.8).  Different 
markers  correspond  to  different  geophones  (see  legend).  The  best  least-square  fit  lines 
for  different  geophones  are  shown  by  color-coded  (the  color  of  the  line  is  the  same 
as  a  color  of  the  marker  corresponding  to  the  respective  geophone)  dashed  lines.  For 
correspondence  between  geophone  names  and  geophone  numbers  used  throughout  the 
text  see  Fig.  3-8. 


Meters  North 


Figure  3-8:  Geophone  array  layout.  The  geophones  are  shown  by  filled  triangles  with 
geophone  names  written  near  them.  Numbers  in  parenthesises  are  used  to  reference 
geophones  throughout  the  text.  So,  for  example,  geophone  G20  is  referred  to  as 
geophone  #3  in  the  text. 


Figure  3-9;  Processing  results  of  the  Event  F5-5  (see  in  the  text)  with  the  HV  MPD 
method.  For  definitions  of  Xmpd  and  Fmpd  see  Equations  (3.5)  and  (3.6).  Different 
markers  correspond  to  the  results  of  processing  of  different  geophones  (see  legend). 
For  correspondence  between  geophone  names  and  geophone  numbers  used  throughout 
the  text  see  Fig.  3-8. 


Figure  3-10:  Processing  results  of  the  Event  F5-5  (see  in  the  text)  with  the  HiV 
MPD  method  are  shown.  For  definitions  of  .X’mpd  and  Impd  see  Equations  (3.7)  and 
(3.8).  Different  markers  correspond  to  the  results  of  processing  of  different  geophones 
(see  legend).  The  best  least-square  fit  lines  for  different  geophones  are  shown  by  color- 
coded  (the  color  of  the  line  is  the  same  as  a  color  of  the  marker  corresponding  to  the 
respective  geophone)  dashed  lines.  For  correspondence  between  geophone  names  and 
geophone  numbers  used  throughout  the  text  see  Fig.  3-8. 
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the  other  (further  referred  to  as  Sg)  -  is  the  internal  Cartesian  coordinate  system  of 
each  geophone^. 

For  the  coordinate  system  §o  the  direction  of  the  x-axis  is  due  East  and  that 
of  the  y-axis  is  due  North  (see  top  plot  on  Fig.  3-5).  For  the  coordinate  system 
§j  the  direction  of  the  x-axis  is  due  North  and  that  of  the  y-axis  is  due  West  (see 
bottom  plot  on  Fig.  3-5).  Hence,  in  order  to  be  able  to  use  the  results  of  the  motion 
product  detector  for  source  location  triangulation,  we  need  to  reconcile  the  different 
conventions  in  these  two  coordinate  systems.  For  triangulation  purposes  it  is  better 
to  choose  the  convention  of  Sa,  so  for  the  HV  MPD  processor,  actually  is 

plotted  versus 

(3.5) 

(3.6) 

where  Xuy  and  Vhv  are  defined  by  Equations  (3.1)  and  (3.2),  respectively. 

Similarly,  for  the  HiV  MPD  processor  Fmpd  is  plotted  versus 

^MPD  =  -^HiV  (3.7) 

=  XhIVj  (3*^) 

where  Xhiv  and  Y^y  are  defined  by  Equations  (3.3)  and  (3.4),  respectively. 


-^MPD  “ 

VHV  _  Y 
^MPD  “ 


3.4  Results  of  data  processing  by  the  MPD  method 

The  processed  data  was  collected  on  a  small  ice  floe  which  was  located  approximately 
4  km  North-East  of  the  main  camp.  The  geophone  array  layout  at  that  site  and 
the  correspondence  between  geophone  names  (used  on  plot  legends)  and  geophone 
numbers  (used  throughout  the  text)  are  shown  on  Fig.  3-8. 

Here  I  will  consider  in  detail  the  results  of  the  two  following  events: 

^  Due  to  the  “conceptual”  rotation  procedure  (see  above)  the  orientations  of  the  main  axes  of  all 
geophones  are  identical  within  the  precision  of  angular  measurement  of  geophones’  orientations 
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1.  This  event  occured  approximately  325  sec.  after  the  beginning  of  file  #5  on  the 
RLAM  tape  #26  (further  referred  to  as  RLAM-26).  This  event  will  be  further 
referred  to  as  Event  F5-3. 

2.  The  second  event  occured  approximately  390  sec.  after  the  beginning  of  RLAM- 
26.  This  event  will  be  further  referred  to  as  Event  F5-5. 

3.4.1  MPD  seismograms 

On  the  Fig.  3-6  the  processing  results  of  the  Event  F5-3  with  the  HV  MPD  method 
are  shown  (recall  that  in  this  processor  there  is  no  phase  shift  between  the  horizontal 
and  vertical  components).  The  different  kinds  of  markers  correspond  to  different 
geophones.  On  this  plot  we  see  points  scattered  all  over  the  place;  there  is  no  obvious 
linear  dependence  between  them.  That  means  that  no  wave  with  strong  rectilinear 
polarization  is  present  in  the  data.  On  the  Fig.  3-7  the  processing  results  for  the 
same  event  with  the  HiV  MPD  method  are  shown  (recall  that  in  this  processor 
there  is  a  phase  shift  of  7r/2  between  the  horizontal  and  vertical  components).  Again, 
the  different  kinds  of  markers  correspond  to  the  results  for  different  geophones.  The 
results  are  strikingly  different  from  that  presented  on  Fig.  3-6.  Really,  we  have  quite 
evident  linear  dependence  in  the  output  for  each  geophone.  The  best  least-square  fit 
lines  for  different  geophones  are  shown  by  color-coded  (the  color  of  the  line  is  the 
same  as  a  color  of  the  marker  corresponding  to  the  respective  geophone)  dashed  lines. 
These  best  fit  lines  are  also  shown  on  separate  plots  for  each  geophone:  see  Figs.  3-11, 
3-12,  3-13,  3-14,  and  3-15.  The  correlation  coefficients  for  those  fits  are  quite  high 
(usually,  above  90%).  Note,  that  lowest  correlations  correspond  to  two  geophones  (#2 
(G9)  and  #3  (G20))  that  give  quite  contradictory  localization  results  when  compared 
to  other  three  geophones,  although  the  actual  difference  in  correlation  coefficients  is 
not  really  significant  (for  example,  94.3%  -  for  geophone  #3  (G20)  and  94.8%  -  for 
geophone  #4  (GIO)). 

A  similar  conclusion  could  be  reached  by  looking  at  Figs.  3-9  and  3-10  which 
show  processing  results  for  the  Event  F5-5  by  the  HV  MPD  and  HiV  MPD  meth- 
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ods  respectively.  Note  that  these  two  events  are  separated  by  at  least  1  min  in  time, 
but  still  the  similarities  between  them  are  quite  evident. 

3.4.2  Triangulation  results 

Fig.  3-16  shows  the  best  fit  lines  drawn  from  corresponding  geophone  locations  for  the 
Event  F5-3  processed  by  the  HiV  MPD  method.  On  this  plot  the  source  location 
could  be  found  by  successful  triangulation  with  at  least  3  bearings  from  geophones 
(the  result  of  such  triangulation  is  shown  by  filled  circle  on  Fig.  3-16).  For  the  other 
two  geophones  (note  that  for  these  two  geophones  linear  fits  show  lower  correlation 
coefficients  than  that  for  other  three  geophones)  bearings  are  quite  contradictory  to 
the  triangulated  source  location.  A  similar  plot  for  the  Event  F5-5  (see  Fig.  3-17) 
shows  identical  behavior:  again  geophones  #1  (GO),  #4  (GlO)  and  #5  (G14)  give 
good  triangulation  for  the  source,  while  geophones  #2  (G9)  and  #3  (G20)  contradict 
them.  A  possible  explanation  for  the  difference  in  the  behavior  of  these  two  geo¬ 
phone  groups  is  the  hypothesis  that  the  ice  plate  could  have  been  broken  somewhere 
between  these  two  groups  and  as  a  result  geophones  #2  and  #3  might  have  been 
on  a  different  ice  plate  than  the  other  three  geophones.  The  horizontal  refraction 
on  the  hypothetical  boundary  of  two  plates  could  account  for  differences  in  apparent 
locations  of  ice  events  for  these  two  geophone  groups. 

The  above  procedure  has  been  applied  to  several  ice  events.  One  result  of  this 
polarization  analysis  was  that  these  fractures  seemed  to  generate  vertically  polarized 
shear  (SV)  waves  for  the  most  part.  The  source  distances  from  the  array  center  for 
these  events  are  shown  versus  time  on  Fig.  3-18.  The  analysis  of  this  figure  indicates 
that  stresses  redistribute  and  move  away  from  the  corner  of  two  ridges  with  time. 

On  Fig.  3-19  the  same  events  are  shown  overlayed  on  the  aerial  photograph  of 
the  experimental  site.  On  this  figure  we  can  see  that  events  mostly  concentrate  in  the 
vicinity  of  the  corner  of  two  ice  ridges.  This  indicates  the  area  of  high  stress.  Also  the 
character  of  particle  motion  in  the  corresponding  arrivals  supplies  additional  evidence 
for  the  continuing  ridge  building  process  in  that  part  of  the  ice  island,  because  ice 
breaking  during  ridge  building  is  characterized  by  predominantly  vertical  particle 
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Figure  3-16:  Results  of  the  triangulation  of  the  source  location  by  best  least-square 
fit  lines  for  the  Event  F5-S  processed  by  the  HiV  MPD  method.  The  triangulated 
source  location  is  shown  by  a  filled  circle.  The  geophones  #1  (GO),  #4  (GlO)  and  #5 
(G14)  produce  consistent  results,  while  the  geophones  #2  (G9)  and  #3  (G20)  give 
quite  different  results. 
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Figure  3-17:  The  results  of  the  triangulation  of  the  source  location  by  best  least-square 
fit  lines  for  the  Event  F5-5  processed  by  the  HiV  MPD  method.  The  triangulated 
source  location  is  shown  by  a  filled  circle.  Again,  like  for  the  Event  F5-3  the  geophones 
#1  (GO),  #4  (GIO)  and  ^5  (G14)  produce  consistent  results,  while  the  geophones 
#2  (G9)  and  #3  (G20)  give  quite  different  results  (compare  with  Fig.  3-16). 
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Figure  3-18:  Ranges  from  array  apex  versus  time  for  events  located  by  triangulation 
using  HiV  MPD  processing  results. 


motion  as  opposed  to  the  horizontal  particle  motion  in  compressional  waves.  The 
stress  relief  is  most  probably  occurring  in  the  ice  lead  which  could  be  identified  on 
the  figure  by  the  location  of  hydrophone  deployment. 


3.5  Conclusions 

Polarization  methods  are  well  known  in  seismology,  but  they  have  never  been  used  for 
ice  event  data  processing.  In  this  chapter  one  of  the  polarization  methods  (so  called 
Motion  Product  Detector  method)  has  been  successfully  applied  for  the  localization 
of  ice  events  recorded  by  the  REAM  box  deployed  on  “ice  island” ,  which  is  a  remote 
experimental  site,  approximately  4  km  East  of  the  base  camp,  and  the  determination 
of  the  polarization  characteristics  of  the  elastic  waves  generated  by  these  events.  The 
fractures  identified  by  this  method  seemed  to  generate  vertically  polarized  shear  (SV) 
waves  for  the  most  part.  Their  concentration  in  the  vicinity  of  the  corner  of  two  ice 
ridges  and  the  character  of  particle  motion  in  the  corresponding  arrivals  indicate  the 
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•  Hydrophone  ♦  Event  location 

Figure  3-19:  Event  locations  determined  by  triangulation  using  HiV  MPD  pro¬ 
cessing  results  overlayed  on  the  aerial  photo  of  the  “ice  island”.  On  this  figure  the 
geophone  locations  are  shown  by  filled  triangles  with  geophone  names  and  reference 
numbers  written  next  to  markers.  The  location  of  hydrophone  is  shown  by  filled 
circle. 
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continuing  ridge  building  process  in  that  part  of  the  ice  island. 

This  application  demonstrated  the  feasibility  of  the  polarization  method  for  ice 
event  data  processing,  because  it  allows  to  identify  areas  of  high  stress  concentration 
and  “hot  spots”  in  ridge  building  process.  Furthermore,  the  polarization  characteris¬ 
tics  of  elastic  waves  could  be  used  for  establishing  the  physics  of  the  processes  in  the 
ice  sheet  which  lead  to  generation  of  ice  events. 
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Chapter  4 


Edge  waves 

4.1  Overview 

In  several  ice  events  recorded  on  the  so-called  “ice  river”  site,  the  site  near  the  ice  lead 
approximately  2  km  north  of  the  base  camp,  evidence  was  found  for  the  existence  of 
so-called  ‘edge  waves’:  waves  propagating  along  the  edges  of  a  newly  opened  lead. 

The  existence  of  similar  ‘edge  wave’  phenomena  have  been  reported  in  the 
literature.  Recently,  Goldstein  et  al.  [14]  demonstrated  theoretically  the  existence  of 
a  new  type  of  waves  in  the  water/ice  system  containing  a  rectilinear  crack  extending 
through  the  ice  cover.  They  found  that  edge  waves  were  propagating  in  the  water 
along  the  crack  and  inducing  the  ice  cover  deformation.  The  role  of  the  ice  sheet 
in  their  study  was  to  modify  the  boundary  conditions  at  the  upper  boundary  of  the 
fluid,  i.e.  the  ice  was  only  involved  as  an  external  loading.  This  loading  (mainly 
stiffness,  because  they  considered  an  inertia-less  thin  plate  model  of  the  ice  sheet) 
modified  the  dispersion  relation  for  the  gravity  surface  waves  on  the  sea  surface.  The 
crack  considered  was  essentially  a  line  of  contact  between  two  ice  plates  which  were  not 
interacting  in  bending.  Goldstein’s  solution  describes  the  diffraction  of  the  waterborn 
bending-gravitational  waves  at  the  infinitely  thin  ice  crack.  However,  Goldstein’s 
approach  does  not  directly  apply  to  the  SIMI-94  experimental  observations:  first,  by 
using  the  model  for  the  ice  sheet  as  an  inertia-less  thin  plate  it  discards  all  types  of 
elastic  waves  that  propagate  in  the  ice  sheet;  further,  the  model  of  the  crack  as  having 
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zero  width  certainly  can  not  be  applied  to  the  ice  lead. 

An  extensive  literature  exists  on  waves  in  plates  without  fluid  loading,  with 
some  papers  devoted  to  ‘edge  waves’  phenomena.  Shaw  [47]  found  experimental  evi¬ 
dence  for  ‘edge  waves’  in  thick  barium  titanate  disks.  He  investigated  their  resonant 
properties  with  the  objective  of  using  them  in  transducers.  He  discovered  that,  as 
the  frequency  approached  a  certain  value,  depending  on  the  disk  thickness,  the  disk 
vibration  was  consistent  with  a  surface  wave  resonance  with  maximum  motion  occur¬ 
ring  at  the  edge  of  the  disk.  He  also  found  that  this  resonance  was  formed  by  several 
modes  of  the  elastic  disk. 

An  attempt  to  describe  this  experimental  data  theoretically  was  first  made  by 
Gazis  and  Mindlin  [13].  They  set  up  approximate  equations  of  extensional  motion  of 
an  elastic  plate  using  a  modal  expansion  in  terms  of  a  finite  number  of  the  modes 
of  an  infinite  plate.  Their  theory  was  limited  to  axially  symmetric  vibrations  of  a 
circular  disk.  They  showed  that  the  ‘edge  mode’  in  this  case  arises  from  the  complex 
conjugate  roots  of  the  frequency  equation  for  the  modes  of  the  infinite  plate.  They 
investigated  the  properties  of  the  ‘edge  mode’  in  the  study  of  the  reflection  of  straight- 
crest  extensional  waves  at  the  edge  of  a  semi-infinite  plate.  In  this  geometry  the 
‘edge  mode’  was  revealed  by  the  substantially  increased  amplitudes  (up  to  7.7  times 
the  amplitude  of  the  incident  wave)  of  the  reflected  complex  waves.  Nevertheless, 
Gazis  and  Mindlin  did  not  achieve  quantitative  agreement  with  the  experimental 
data  obtained  by  Shaw.  As  shown  later  by  Torvik  [51],  the  main  reason  was  that 
Gazis  and  Mindlin  used  an  approximation  of  the  frequency  spectrum  rather  than 
exact  roots  of  the  Rayleigh-Lamb  equation,  and  too  few  modes  were  employed  in  the 
model.  Torvik  improved  the  model  by  using  the  exact  frequency  spectrum  and  by 
using  many  more  modes,  and  a  greatly  improved  agreement  with  Shaw’s  experiment 
was  achieved. 

The  more  recent  work  by  Gregory  and  Gladwell  [15]  used  a  ‘method  of  pro¬ 
jections’  developed  by  the  authors  for  solving  similar  problems  in  elastostatics.  They 
represented  the  stress  function  as  a  sum  of  some  inhomogeneous  term  and  series  of 
orthogonal  functions.  Their  approximate  solution  was  the  best  approximation  in  the 
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least-square  sense,  and  they  again  confirmed  the  existence  of  the  ‘edge  mode’  at  a 
single  isolated  value  of  the  frequency  by  considering  the  normal  incidence  of  the  plain- 
strain  symmetric  Rayleigh-Lamb  wave  onto  the  free  edge  of  the  semi-infinite  plate. 
They  found  that  the  resonance  peak  was  much  narrower  and  higher  than  previously 
found  by  Torvik,  and  Gazis  and  Mindlin.  However,  the  peak  was  bounded,  indicating 
energy  dissipation  in  the  near-edge  region. 

A  common  feature  of  all  of  this  past  work  is  the  resonant  nature  of  the  ‘edge 
mode’,  appearing  at  discrete  frequencies.  Also,  the  past  theoretical  work  considered 
normal  incidence  waves  on  the  free  edge.  This  is  significantly  different  from  what  was 
observed  during  SIMI-94:  the  waves  observed  along  the  lead  were  of  true  broadband 
nature,  and  similar  in  behavior  to  normal  seismic  waves  in  the  ice  sheet,  but  with  a 
characteristic  periodicity  and  propagation  direction  parallel  to  the  lead.  The  reason 
for  the  existence  of  “edge  mode”  is  that  a  flexural  wave  propagating  along  a  stress- 
release  ice  edge  will  be  slower  (the  plate  is  softer  at  the  edge)  than  a  flexural  wave  in 
the  infinite  plate.  This  speed  gradient  near  the  edge  provides  a  waveguide. 

The  problem  in  the  theoretical  treatment  of  the  problem  is  the  free  edge  of 
the  ice  plate  which  makes  the  problem  non-separable.  Therefore,  in  contrast  to  the 
infinite  plate  problem  [41],  [39],  a  simple  spectral  integral  solution  cannot  be  derived. 

These  circumstances  led  to  the  construction  of  an  approximate  model  which 
could  explain  the  characteristics  and  behavior  of  the  observed  edge  waves.  The  model 
was  based  on  the  lowest  order  modes  of  the  infinite  plate,  coupled  through  integral 
boundary  conditions  at  the  free  edge. 

This  chapter  first  describes  the  experimental  setup,  followed  by  the  process¬ 
ing  schemes  applied.  Then  the  approximate  theoretical  model  is  described,  and  the 
comparison  with  the  experimental  data  is  discussed. 
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Figure  4-1:  Satellite  synthetic  aperture  radar  images  of  SIMI-94  site.  The  upper  plot 
was  recorded  at  7.11  am  GMT  on  April  17,  1994.  The  lower  plot  was  recorded  at 
9.48  pm  GMT  on  April  18, 1994.  The  yellow  dots  near  the  center  of  each  image  show 
the  location  of  the  main  hydrophone  array  (see  Fig.  2-2).  The  new  lead  is  evident  as 
a  darker  NW-SE  strip  to  the  NE  of  the  array  on  the  lower  plot,  while  there  is  nothing 
remarkable  in  the  corresponding  part  of  the  upper  plot. 
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Figure  4-2:  High-resolution  satellite  synthetic  aperture  radar  image  of  SIMI-94  site 
recorded  at  9.48  pm  GMT  on  April  18,  1994.  The  black  dots  near  the  center  of  the 
image  show  the  location  of  the  main  hydrophone  array  (see  Fig.  2-2).  The  new  lead 
is  evident  as  a  darker  NW-SE  strip  to  the  NE  of  the  array,  with  the  local  direction 
indicated  by  the  solid  line.  Letters  A  and  B  indicate  the  approximate  locations  of 
the  geophone  arrays. 
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4.2  The  experimental  setup  and  the  geophone  data 

The  data  in  which  evidence  for  edge  waves  was  found  was  collected  in  the  vicinity 
of  a  major  lead  formed  April  17-18,  1994.  This  can  be  seen  from  two  SAR^  satellite 
images  from  the  Alaskan  SAR  Facility  taken  on  those  two  dates  (see  the  Fig.  4-1). 
This  lead  can  be  observed  in  the  upper  right  corner  of  the  lower  plot  of  this  figure 
(taken  on  April  18),  while  there  is  nothing  remarkable  in  the  corresponding  part  of 
the  upper  plot.  The  closer  look  of  this  lead  is  given  on  Fig.  4-2  (also  recorded  on 
April  18).  On  this  figure  the  predominantly  SE-NW  direction  of  the  lead  is  indicated 
by  a  solid  line,  and  the  surveillance  array  position  is  superimposed  as  well.  At  its 
closest  point,  the  distance  from  the  camp  to  the  lead  is  approximately  2  km.  Two 
REAM  geophone  arrays,  each  having  5  geophones  and  1  hydrophone,  were  deployed 
in  the  immediate  vicinity  of  the  lead  on  its  SW  edge,  as  indicated  in  Fig.  4-3.  Each 
geophone  measured  three  components  of  velocity  of  ice  surface  movement,  and  the 
hydrophone  measured  the  acoustic  pressure  in  the  water.  The  aperture  (radius)  of 
the  two  geophone  clusters  was  approximately  50  m,  and  the  separation  along  the  lead 
was  400  m. 

Fig.  4-4  shows  data  from  2  geophones  in  the  site  A  cluster.  The  quasi¬ 
periodicity  of  the  signal  is  evident  both  on  the  time  series  plot  and  on  the  autocorre- 

^  SAR  stands  for  Synthetic  Aperture  Radar 
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Figure  4-3:  Geophone  arrays  A  and  B  deployed  on  SW  shore  of  new  lead  2  km  NE  of 
camp.  Hydrophones  were  deployed  in  the  lead  at  a  depth  of  60  m.  The  approximate 
directions  of  the  3  geophone  axes  are  indicated  by  the  coordinate  system  {X,Y,Z). 
The  aperture  (radius)  of  both  geophone  clusters  was  approximately  50  m,  and  the 
separation  between  them  along  the  lead  was  400  m. 

lation  plot  for  the  Y-component  of  the  ice  particle  velocity  recorded  by  the  geophone 
#lb  from  site  A  (see  Fig.  4-5).  The  frequency- wavenumber  processing  described 
below  was  performed  event  by  event,  with  the  event  intervals  specified  in  Table  4.1. 

4.3  Data  processing  scheme 

The  wavenumber  processing  is  based  on  conventional  plane  wave  beamforming  using 
estimated  covariance  matrices.  The  sample  covariance  matrix  Sx(ti^)  at  angular  fre¬ 
quency  U3  is  obtained  from  the  set  of  complex  Fourier  component  of  the  signal,  A|(a;), 
as 

..  M 

/=! 

where  z,  j  represent  the  sensor  numbers,  I  is  the  segment  number  (/  <  M)  and  the 
asterisk  denotes  complex  conjugation.  To  obtain  each  frequency  component,  each 
segment  of  the  raw  data  is  tapered  by  a  Hamming  window  and  Fourier  transformed 
using  an  FFT  algorithm.  To  avoid  energy  leakage  at  the  edges  of  the  taper  window, 
the  segments  are  chosen  to  overlap  by  50%.  Now  the  estimate  of  the  power  spectrum 
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Figure  4-4:  Outputs  of  2  geophones  from  site  A.  The  labeling  (X,Y,Z)  represents  the 
geophone  components  consistent  with  Fig.  4-3  (X-axis  is  towards  the  lead,  Y-axis  - 
parallel  to  the  lead).  Event  intervals  are  specified  in  Table  4.1. 


93 


-Too 


-50 


0 

Lag,  sec. 


50 


100 


Figure  4-5:  Autocorrelation  for  the  Y-component  of  the  ice  particle  velocity  recorded 
by  the  geophone  #lb  from  site  A.  The  time  interval  (from  40  sec.  to  120  sec.)  covers 
all  8  events  from  table  4.1. 
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at  frequency  ui  is 


PconvCw,  k)  =  (^E+  /N^  ,  (4.1) 

with  E  being  the  array  steering  vector 

E  =  E{k)  =  : 

^  0-jk'^XN 

Here  k  is  the  horizontal  wavenumber  vector  in  the  ‘look’  direction,  f*  is  the  location  of 
the  i-th  sensor,  and  N  is  the  total  number  of  sensors.  Each  component  of  the  geophone 
output  is  processed  separately,  which  implies  N  —  ^  (number  of  geophones  in  each 
array).  For  each  frequency  the  estimated  2-D  wavenumber  spectrum  is  contoured, 
and  the  maxima  corresponding  to  the  various  arrivals  are  identified.  Fig.  4-6  shows 
a  characteristic  example  at  one  of  the  frequencies  (15  Hz)  used  in  the  processing, 
showing  two  predominant  propagation  directions,  towards  the  SE  and  towards  the 
NW,  consistent  with  the  lead  orientation. 

Particularly  strong  evidence  for  the  ‘edge  wave’  phenomenon  was  found  in 
events  #2  (Fig.  4-6)  and  #3  at  site  A.  The  frequency  dependence  of  the  propagation 
direction  for  the  two  events  are  shown  in  Fig.  4-7.  The  directions  are  obviously 
identical  to  the  propagation  directions  back  and  forth  along  the  lead,  indicated  by 
the  solid  lines.  The  following  numbers  confirm  that  conclusion.  For  the  event  #3 
(crosses)  the  mean  direction  is  m3  =  126.5°  with  =  5.5°^,  while  for  the  event  #2 
-  the  mean  m2  =  -47.5°  and  ^2  =  2.5°.  We  see  that 

|180°  -  (m3  -  m2)|  <  (<73  -I-  (T2). 


^  (7  is  the  square  root  of  standard  deviation  about  the  mean  value 
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Figure  4-7:  Directions  of  wave  arrivals  for  ice  events  #2  and  #3.  The  y-axis  indicates 
the  angle  a  of  propagation,  counter-clockwise  from  due  East.  Dashed  lines  mark 
the  mean  angle  of  propagation,  while  the  solid  lines  mark  the  prevailing  local  lead 
direction  as  indicated  in  Fig.  4-2. 
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4.4  Theoretical  model 


4.4.1  General  equations  of  motion  in  elastic  medium 

The  properties  of  the  operator  V 

We  will  need  most  of  the  properties  of  the  operator  V  for  the  following  derivation: 


V  -  (v?i)  =  vV 

(4.2) 

V(V  ■^)  =  VH  +  V  x^) 

(4.3) 

V  X  {V4>)  =  0 

(4.4) 

V  •  (V  X  ijf)  =  0 

(4.6) 

The  original  equation  of  motion  in  homogeneous  elastic  medium 
The  original  equation  of  motion  in  homogeneous  elastic  medium  is 

(A  + /x)V(V  •  u)  +  ^V^u  —  pii  =  0  (4.6) 

where 

u  =  displacement  vector, 

p  =  the  density  of  the  medium, 

A,  /i  =  Lame  coefficients  of  the  medium. 

Using  Equation  (4.3)  we  get 

V^u  =  V(V  -  u)  -  V  X  V  X  u 
and  the  new  form  of  the  equation  of  motion  is 

(A  +  2)u)V(V  •  u)  -  //V  xVxu-pu  =  0  (4.7) 
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The  potentials 

The  displacement  field  (as  any  vector  field)  can  be  described  by  a  scalar  potential 
and  a  vector  potential: 

u  =  V(j)  +  V  X  (4-8) 

In  Equation  (4.8)  the  vector  potential  ^  should  satisfy  the  gauge  condition: 

V  •  ^  =  0.  (4.9) 


®  compare  fourth  term  in  Equation  (4.10)  and  second  term  in  Equation  (4.11) 
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pressional  (cp)  and  shear  (cg)  waves 


we  get  the  following  form  of  equation  of  motion 

CpPV(VV  -  4^)  +  CspV  X  -  4'^)  =  0.  (4.13) 

The  solution  to  this  equation  will  be  given  by  the  potentials  0  and  that  satisfy  the 
following  equations: 

VV  -  4*^  =  0  (4.14) 

_  if  =  0.  (4.15) 

For  harmonic  field  with  frequency  u  we  have 


0  =  — 
f  = 


and  Equation  s(4.14)  and  (4.15)  take  the  form  of  Helmholtz  equation: 

0 

0. 

The  vector  displacement  potential  ’5'  can  be  represented  as  follows: 


^  =  V  X  (0,0,  A) +  (0,0,^,). 


(4.16) 


Introducing  compressional  media  wavenumber  (kp)  and  shear  media  wavenumber  (kg) 
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as 


(4.17) 


UJ 

h>  —  _ 

rvg  —  5 


(4.18) 


we  obtain  the  following  form  of  the  Helmholtz  equations  for  the  displacement  poten¬ 
tials  of  the  homogeneous  part  of  the  field 


+  (l){x,y,z)  =  0 

(V^  -h  k^)  Aix,  y,z)^0 

(V^  +  k'^s)  '^z{x,y,z)  =  0. 

Let  us  apply  the  Fourier  Transform  with  respect  to  y 


/+0O 

f{x,  y,  z)e'^^ydy 

■oo 

1  /'+0O 

f{x,  yy^)  =  -^  J 


on  the  Helmholtz  equations  to  get: 


-  '‘D)  =  ^ 

(gj2  5^2  “  ~  *»'  ^  **■ 


(4.19) 

(4.20) 

(4.21) 


(4.22) 

(4.23) 


Similarly,  after  applying  the  Fourier  Transform  with  respect  to  a;,  we  obtain: 


(4.24) 

(4.25) 

(4.26) 


A(k,,k„z)=0 

^-f]-iAk„k,,z)  =  0, 


(4.27) 

(4.28) 

(4.29) 
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where 


«=Y^x  +  ^y-^p  (4.30) 

^  +  -  kl-  (4.31) 

Solving  above  equations  for  the  displacement  potentials  we  get: 

0(^x,  ky,  z)  =  A~  (kx,  fcj,)e“®^  +  A'^(kx,  ky)e°‘^  (4.32) 

^.(A:x,  ky,  z)  =  B-{kx,  ky)e-P^  +  B+(A:^,  ky)e^^  (4.33) 

K{kx,  ky,  z)  =  C-{kx,  ky)e-^^  +  C-^{kx,  ky)e^^  (4.34) 


Expressions  for  displacement  components 
Now  for  displacement  we  have  the  following  expression 


u{x,  y,  z)  =  V(f>{x,  y,  z)  +  V  X  (V  X  (0, 0,  A)  +  (0, 0,  "J^^)) 
=  V^  +  Vx(|-A.-Aa.^.). 

Introducing  the  notation  u  =  Ux,  v  =  Uy,  w  =  and 

=  A-{kx,  ky)e-^^  ±  A+(kx,  ky)e^^ 

=  B-{kx,  ky)e-P^  ±  B-^{kx,  ky)e^^ 

=  C-{kx,  ky)e-^^  ±  C+{kx,  ky)e^^ 

we  get 

dcj)  d^A  .  u.  .  ^  . 

dx'^'d^'^  ~d^z  ~  ~  +  ^kxPSc 

d(l)  d^A  .  a.  .  ^ 

dy~  ~  +  ikyPSc 

d(f)  /  d^\  ^  ^ 


(4.35) 


(4.36) 

(4.37) 

(4.38) 
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Expressions  for  stresses 


Using 


and 


we  have  for  stresses: 


/X  r.  \ 

(4.39) 

^ZX  —  M  j 

f  dw 

\dx  ^  dz) 

(4.40) 

II 

( dw  dv\ 

\dy  ^  dz) 

(4.41) 

/X  rw  A  ( d'w\ 

axx  -  (A  -t-  2/x)—  +  + 

(4.42) 

^xy  —  j 

fdu  \ 

{dy  dx) 

(4.43) 

=  A*  1 

( du  dw\ 

\dz  ^  dx  ) 

(4.44) 

Note,  that  of  course,  is  always  equal  to  o^z-  Both  of  them  are  kept  in  the 
derivation  for  convenience  of  formulating  the  boundary  condition  at  the  vertical  free 
face  of  the  plate  at  a:  =  0. 

Substituting  Equations  (4.36-4.38)  into  above  expressions  and  introducing  the 
following  notation  for  horizontal  wavenumber 
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we  obtain: 


=  (A  +  2//)  ^  (_^25+  _  ^ 

+  A  {-klSt  +  KkyS+  +  klpSc) 

=  /i(27^  -  A:^)5+  -  2n'y'^pSc, 

Ozx  =  A*  {ikxoS^  -  ikx'fS'^)  +  /i  -  ik^^^S'^) 

=  2i^ikxQS^  +  inkyPSg  -  inkx{2'f  -  k^)S^, 
azy  =  11  (ikyaS^  -  ikyj^S+)  +  fi  {ikyaS^  -  ikyfS^  -  i^^Sg) 

=  2iiikyaS^  -  i/ikxPSg  -  ifiky{2'y^  -  k‘^)S^, 

(^xx  =  (A  +  2//)  (^—kj.S^  —  kxkySg  +  + 

+  ^  +  kxkyS+  +  klpSc)  +  A  {a^St  - 

=  [-*x(^  +  2a^)  -  AA:2  +  Aa^]  S+  +  [-k^kyiX  +  2fi)  +  Xk^ky]  S+  + 

+  [A:|^(A  +  2/i)  +  Xk^P  -  Xj^p]  Sc 

=  -{2fikl  +  AA:^)5^  -  2(ikxkyS+  +  2tiklpSc, 

Oxy  =  fi  {-kxkyS+  -  klS+  +  kxkypSc)  +  n  {-kxkyS+  +  klS+  +  kxkypSc) 
=  —2fikxkyS'^  +  ti{k^  —  ky)Sg  +  2^kxkyPSc, 

Oxz  —  ifikxCxSj^  +  ifikyPSg  ifikxP^ Sq  +  fi  (ikxCxSy^  —  ikxy^ Sq') 

=  2ifikxaS2  +  i^ikypSg  -  ifikx{2'y^  -  k1)S^. 
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Now  we  have  the  final  expressions  for  stresses: 


Ozz  =  (27^  -  (4.45) 

(72*  =  2iiika;aS2  +  ifjikyPSg  -  i^k^  (27^  -  kfj  5^  (4.46) 

Ozy  -  2ifj,kyaS^  -  ifik^PS^  -  ifxky  (27^  -  A:^)  (4-47) 

a**  =  -  (2^kl  +  \kl)  S\  -  2fik,kyS+  +  2fiejSc  (4.48) 

(7*y  =  -2fikxkyS'^  +  n{kl  -  kl)S^  +  2nk^kypSc  (4.49) 

a*2  =  2ii^kxaS2  +  i/iky^Sg  -  inkx{2'f  -  kl)S'^  (4.50) 


We  will  need  Eqs.4.36-4.38  and  Eqs.4.45-4.50  to  write  expressions  for  boundary  con¬ 
ditions  at  the  free  vertical  face  of  the  plate. 


4.4.2  The  fundamental  modes  of  the  infinite  plate 

As  a  basis  for  the  edge  wave  model,  the  first  symmetric  and  antisymmetric  Lamb 
modes  (P-SV)  and  first  symmetric  and  antisymmetric  Love  modes  (SH)  of  the  infinite 
plate  are  assumed  to  be  dominating.  This  is  a  reasonable  assumption  for  frequencies 
below  the  t/iicfcness-sftear  frequency  for  the  the  ice  cover,  typically  of  the  order  200-300 
Hz.  The  discussion  of  these  modes  follows  below. 

The  first  symmetric  Lamb  mode 

For  the  first  symmetric  Lamb  mode  we  have  the  following  symmetry  conditions; 

A+  =  A-  ^  A 
B+  =  B-  =  0 
(7+  =  -C~  ^  C. 
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These  give  the  following  depth  dependencies 


=  2Aco8\i{az) 

=  — 2Asinh(Q:^) 

5+  =  0 

5^  =  0 

S'^  =  2Csm\i{^z) 

= —2C  cosh{Pz) 

In  the  above  expressions  constants  A  and  C  could  be  expressed  in  terms  of  one 
constant  i4syin  using  one  of  the  boundary  conditions  on  either  one  of  the  horizontal 
faces  of  the  plate  (either  aX  z  =  -h  ox  zX  z  =  h): 

A  =  27^^  cosh(y0h)zlsyni 
C  =  -(27^  -  k1)  cosh(o;/i)Asyni. 

The  above  depth  dependencies  give  us  following  expressions  for  the  stress 
components  which  are  used  for  the  vertical  boundary  condition  at  x  =  0: 

(7^^  =  -2  {2nkl  +  AA:J)  A  cosh(«2)  -  A^kl^C  cos\i{^z)  (4.51) 

(^xy  =  -^jxkj^kyA  cosh (02:)  -  Ankj^kyPC cosh{Pz)  (4.52) 

(^xz  =  -4iij,kj;aAsmh{az)  -  2ifik3;{2j^  -  kg)C smh{l3z)  (4.53) 

The  dispersion  relation  for  the  first  symmetric  Lamb  mode  could  be  easily  obtained 
from  the  symmetry  conditions  stated  above  and  the  boundary  conditions  at  either 
one  of  the  horizontal  faces  of  the  place  (either  sX  z  =  — h  or  sX  z  =  K).  The  derivation 
of  the  dispersion  relation  will  be  given  only  for  this  mode  (for  other  modes  it  could 
be  done  by  analogy).  We  have  two  equations  zX  z  =  h: 
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This  gives  us  the  following  system  of  equations: 


2//(27^  —  k^)Acosh{ah)  +  I3C  cosh{Ph)  =  0, 

— 4i/xA:a;Q;^sinh(Q!/i)  —  2i^kx{2'y^  —  kl)C  smh{Ph)  =  0. 

After  dividing  the  first  one  by  2/i  and  the  second  one  by  —2ifikx  we  get 

(2j^  —  k^)  cosh{Q;h)A  +  27Vcosh(/0h)  (7  =  0, 
2Q:sinh(Q:h)A  —  (27^  —  k^)  sinh(;5h)(7  =  0. 


This  system  has  nontrivial  solutions  with  respect  to  A  and  C  only  when  the  deter¬ 
minant  of  the  system  is  zero: 


det 


(27^  —  A:|)  cosh  (ah) 
2a  sinh(ah) 


27^/0  cosh  (/3h) 

-(27^  -  k'^^)smh{ph) 


=  0 


(4.54) 


This  equation  leads  to  the  following  form  of  the  dispersion  relation  for  the  symmetric 
Lamb  modes 


tanhySh  _ 

tanh  ah  (27^  —  k^)"^ 


(4.55) 


We  could  solve  it  numerically  with  respect  to  7.  As  a  result,  we  would  obtain  7  as  a 
function  of  frequency  w  after  choosing  branch  corresponding  to  the  first  mode: 


7  =  Vs{uj). 


(4.56) 


The  first  antisymmetric  Lamb  mode 

For  the  first  antisymmetric  Lamb  mode  we  have  the  following  symmetry  conditions: 


= 

C+  = 


-A-^A 
B-  =  0 
c-  =  C. 


These  give  the  following  depth  dependencies 


=  2Asinh(Q;z) 

5^^  =  — 2Acosh(Q!2:) 

<5^  =  0 

5^=0 

Sq  =  2C  cosh(^2;) 

Sq  =  —2(7  sinh(/52:) 

Again,  in  the  above  expressions  constants  A  and  C  could  be  expressed  in  terms  of  one 
constant  Aasym  using  one  of  the  boundary  conditions  on  either  one  of  the  horizontal 
faces  of  the  plate  (either  a-t  z  =  —h  or  at  z  =  h): 

A  =  27^/5  sinh{ph)A^yj^ 

C  =  -(27^  -  k1)  sinh(Q!/i)Aasym. 

The  above  depth  dependencies  give  us  following  expressions  for  the  stress  components 
which  are  used  for  the  vertical  boundary  condition  at  a:  =  0: 

(^xx  =  -2  (2nkl  +  Xkl)  A  sinh(Q;2)  -  Afxkl^C smh{/3z)  (4.57) 

(^xy  =  -‘ifj.ka^ky A  smh{az)  -  Afik^ky/SC  smh{pz)  (4.58) 

CTxz  =  -4ifj,ka;aA  cosh(Q!2)  -  2iixkx{2'f  -  kl)2C  cos\i{pz)  (4.59) 
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The  dispersion  relation  for  the  antisymmetric  Lamb  modes  could  be  easily  obtained 
from  the  symmetry  conditions  stated  above  and  the  boundary  conditions  at  either 
one  of  the  horizontal  faces  of  the  place  similarly,  as  above  for  the  symmetric  Lamb 
modes.  It  has  the  following  form: 

47  V 

tanhah  \{2j^  -  )  ' 

Solving  it  numerically  with  respect  to  7  and  choosing  the  branch  corresponding  to 
the  first  mode  we  get  7  as  a  function  of  frequency  u: 


7  =  Va{u). 


(4.61) 


The  plots  of  the  dispersion  relations  for  the  first  symmetric  and  antisymmetric 
Lamb  modes  of  the  infinite  plate  are  shown  in  non-dimensional  form  on  the  Fig.  4- 
8  where  non-dimensional  frequency  Q  is  plotted  versus  non-dimensional  horizontal 
wavenumber  F.  Non-dimensional  frequency  and  horizontal  wavenumber  F  are  de¬ 
fined  by 


F  = 


2h 

—7 

TT 


The  first  symmetric  Love  mode 

For  the  first  symmetric  Love  mode  we  have  the  following  symmetry  conditions: 

A+  =A-  =  0 
B+  =B-  ^  B 

C+  =c-  =  0. 
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First  two  fundamental  Lamb  modes 


Figure  4-8:  Dispersion  curves  of  the  first  two  fundamental  Lamb  modes  for  Cp  = 
3500  m/s,  Cs  =  1800  m/s  given  in  non-dimensional  form  as  non-dimensional  frequency 
D  versus  non-dimensional  horizontal  wavenumber  F  (two  solid  lines).  Here  Cp  — 
dispersion  for  compressional  wave,  =  dispersion  for  shear  wave,  dashed  line  = 
dispersion  for  the  Rayleigh  wave. 
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These  give  the  following  depth  dependencies 


II 

o 

5;  =  o 

Si=2B 

Sb  =  0 

SJ=0 

O 

11 

The  above  depth  dependencies  give  us  the  following  expressions  for  the  stress  com¬ 
ponents  which  are  used  for  the  vertical  boundary  condition  at  a:  =  0: 


^XX  —  AifXkxkyB 

(4.62) 

(4.63) 

(Txz  =  0 

(4.64) 

So,  the  mode  shape  for  the  first  symmetric  Love  mode  is  a  constant.  The 
dispersion  relation  takes  the  form 


7  —  <5sh(‘*^)  —  ks- 


(4.65) 


The  first  antisymmetric  Love  mode 

For  the  first  symmetric  Love  mode  we  have  the  following  symmetry  conditions: 

A+  =  A-  =  0 
=  -B-  =  B 

C+  =  c-  =  0. 


Ill 


These  give  the  following  depth  dependencies 


0 

II 

Sx=0 

Sg  =  2Bs\nh{pz) 

Sg  =  -2B  cosh(^2:) 

0 

II 

0  1 

II 

0 

The  above  depth  dependencies  give  us  following  expressions  for  the  stress  components 
which  are  used  for  the  vertical  boundary  condition  at  x  =  0: 

<^xx  = -4/ifca;fcj,Bsinh{^2:)  (4.66) 

Cfxy  =  2)U(A:^  -  kl)B  smh{Pz)  (4.67) 

<Jxz  =  —2iij,ky/3B  cosh{Pz)  (4.68) 

For  the  first  antisymmetric  Love  mode  the  dispersion  relation  takes  the  form 

(4.69) 

We  could  solve  it  explicitly  with  respect  to  7.  As  a  result,  we  would  obtain  7  as  a 
function  of  frequency  w: 

7  =  ^SH(a;)  =  (4.70) 

4.4.3  The  geometry  of  the  problem 

The  model  geometry  assumes  a  homogeneous  elastic  plate  with  two  horizontal  faces 
at  2  =  ±h,  and  a  vertical  edge  at  x  =  0  (Fig.  4-9).  This  geometry  is  obviously  not 
separable,  and  does  not  allow  for  closed  form  solutions.  However,  several  approximate 
solutions  to  similar  plate  vibration  problems  have  been  presented. 
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-h, za+h 


Figure  4-9:  Edge  wave  model  geometry.  The  ice  plate  is  semi-infinite  with  horizontal 
faces  at  2:  =  ±/i  and  vertical  edge  at  a:  =  0.  There  is  no  edge  at  y  =  0. 


4.4.4  The  boundary  conditions  in  the  weak  form 

Mindlin  [37]  converts  the  3-D  equations  of  elasticity  to  an  infinite  series  of  2-D  equa¬ 
tions  by  expanding  the  displacement  in  a  power  series  of  the  thickness  coordinate, 
and  integrating  through  the  plate  thickness.  This  depth  integration  is  conceptionally 
similar  to  the  Galerkin  approach  used  here.  In  the  Galerkin  approach  [11],  continuity 
of  the  field  across  the  vertical  boundary  is  expressed  in  the  weak  form: 


a^rdz  =  0 


axydz  =  0. 


(4.71) 

(4.72) 

(4.73) 


Basically,  in  this  approach  we  minimize  the  least-square  error  in  the  vertical 
boundary  condition  across  the  plate  thickness. 

As  a  basis  for  the  edge  wave  model  the  first  symmetric  and  antisymmetric 
Lamb  modes  (P-SV)  and  first  symmetric  and  antisymmetric  Love  modes  (SH)  of  the 
infinite  plate  are  assumed  to  be  dominating.  This  is  a  reasonable  assumption  for 
frequencies  below  the  thickness-shear  frequency  for  the  the  ice  cover,  typically  of  the 
order  200-300  Hz. 

Here,  we  will  consider  the  cases  of  different  modes  incident  on  the  vertical 
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interface  of  the  plate.  These  modes  will  be  the  modes  considered  above.  For  all 
cases  on  the  vertical  interface  we  have  in  all  waves  interacting  at  the  interface  the 
identical  values  of  ky  and  components  of  the  wavenumber  vector  due  to  the  trace 
matching  condition.  We  will  consider  the  y-component  of  the  wavenumber  vector  as 
independent  variable  (denoted  further  as  k)  and  try  to  express  the  other  components 
in  terms  of  k  and  the  medium  wavenumbers.  Note,  the  incidence  is  considered  to  be 
in  the  positive  x-direction,  and  without  fluid  loading  the  plate  geometry  is  symmetric, 
and  the  symmetric  and  antisymmetric  equations  are  uncoupled. 

The  incidence  of  the  first  symmetric  Lamb  mode 

Due  to  the  symmetry  of  the  problem  we  will  have  the  following  waves  reflected  from 
the  vertical  interface  as  a  result  of  the  incidence  on  it  of  the  first  symmetric  Lamb 
mode  (further  referred  to  as  Sj):  reflected  first  symmetric  Lamb  mode  (further  re¬ 
ferred  to  as  Sr),  reflected  first  symmetric  Love  mode  (further  referred  to  as 
The  corresponding  wavenumber  vectors  will  look  like 

Sj  :  k  =  {ki,K,kz) 

Sr  .  k  (  k\,  Ki,  k,,) 

:  k  =  (-fe,  K,  fc.) 


The  horizontal  wavenumber  7  should  be  determined  from  the  corresponding 
dispersion  relation  for  the  first  symmetric  Lamb  mode  (see  Equation  (4.56)): 


7  =  Vs{u). 

We  can  now  determine  the  x-component  ki  of  the  wavenumber  vector  in  the  incident 
wave  from  the  value  of  the  horizontal  wavenumber  and  k  (see  above): 

A:i  =  —  K?. 
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The  vertical  component  of  the  wavenumber  vector  could  be  obtained  as 


kz  =  ^^p-7^- 

After  that  we  can  easily  determine  the  value  of  k2 

k2  =  \/k^  -  {kl  +  K^)  =  \Jk'^s~  +  7^  ~ 

The  incidence  of  the  first  antisymmetric  Lamb  mode 

Due  to  the  symmetry  of  the  problem  we  will  have  the  following  waves  reflected  from 
the  vertical  interface  as  a  result  of  the  incidence  on  it  of  the  first  antisymmetric  Lamb 
mode  (further  referred  to  as  A/):  reflected  first  antisymmetric  Lamb  mode  (further 
referred  to  as  Ar),  reflected  first  antisymmetric  Love  mode  (further  referred  to  as 
SH^yni))  corresponding  wavenumber  vectors  will  be 


A/: 

k  =  {k2,  K,  kz) 

Ar  : 

II 

SH^asym)  . 

II 

The  horizontal  wavenumber  7  should  be  determined  from  the  corresponding 
dispersion  relation  for  the  first  antisymmetric  Lamb  mode  (see  Equation  (4.61)): 


7  =  Va{u}). 


We  can  now  determine  the  ar-component  k2  of  the  wavenumber  vector  in  the  incident 
wave  from  the  value  of  the  horizontal  wavenumber  and  k: 


k2  =  \/7^  - 
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The  vertical  component  of  the  wavenumber  vector  could  be  obtained  as 


K  =  y/kl  -  7". 

After  that  we  can  easily  determine  the  value  of  ki  (though  it  is  not  really  used  in  this 
case) 


*^1  = 

The  solution  of  the  vertical  problem 

Due  to  the  simple  character  of  the  depth-dependence  of  the  lowest  order  modes  (hy¬ 
perbolic  functions)  the  depth  integrals  in  Equations  (4.71-4.73)  could  be  evaluated  in 
the  closed  form.  These  ’weak’  boundary  conditions  then  lead  to  a  system  of  linear 
equations  of  the  form 


GX  =  Y 


where  X  is  the  the  wavefield  amplitude  vector:  X  =  {A,B,C),  Y  represents  the 
contributions  from  an  incident  wave  field. 

The  modal  solutions  are  then  found  from  the  eigenvalue  problem  defined  sim¬ 
ilarly  to  the  homogeneous  field  case  (compare  with  Equation  (4.54)): 

det  G  =  0. 

So  the  procedure  for  the  solution  of  the  vertical  problem  can  be  outlined  as  follows: 

•  Choose  the  incident  mode  from  the  fundamental  Lamb  and  Love  modes  of  the 
infinite  plate. 

•  For  selected  incident  mode  determine  the  components  of  wavenumber  vector  for 
all  waves  interacting  at  the  vertical  interface  in  this  particular  case. 

•  Calculate  contributions  to  axx  and  axy  from  all  the  interacting  waves  after 
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substituting  the  closed  form  expressions  for  the  corresponding  depth  integrals. 

•  Combine  all  the  terms  corresponding  to  reflected  waves  on  the  left  side  of  the  sys¬ 
tem  of  linear  equations  representing  ’weak’  vertical  boundary  conditions  (Equa¬ 
tions  (4.71-4.73))  to  obtain  matrix  G  and  the  reflected  wavefield  vector  X,  while 
on  the  right  side  -  all  terms  corresponding  to  the  incident  mode  to  obtain  the 
vector  Y. 

•  Find  the  modal  solutions  to  the  problem  by  solving  the  equation 

det  G  =  0. 

with  respect  to  y-component  of  the  wavenumber  vector  k. 

4.5  Comparison  of  model  and  the  data 

The  model  results  for  mode  dispersion  have  been  compared  to  the  results  of  the 
experimental  analysis.  The  closest  match  was  obtained  for  the  antisymmetric  prob¬ 
lem,  indicated  by  the  solid  line  in  Fig.  4-10.  Between  10  and  20  Hz  the  theoretical 
prediction  agrees  reasonably  well  with  the  experimental  values.  On  the  other  hand, 
the  frequency  dependence  is  obviously  rather  different  for  the  model  and  the  data, 
suggesting  that  a  more  complete  model  is  required.  The  antisymmetric  modes  have 
significant  vertical  motion  components,  as  is  also  evident  in  the  experimental  data, 
and  it  is  therefore  hypothesized  that  the  fluid  loading  would  significantly  alter  the 
dispersion  relation.  Thus,  the  added  mass  is  expected  to  lower  the  phase  velocity  at 
the  low  frequency  end,  as  required  for  improving  the  agreement  in  Fig.  4-10. 

Additionally,  polarization  processing  discussed  in  previous  Chapter  was  car¬ 
ried  out  for  the  events  involved  in  edge  waves  phenomenon.  The  results  did  not  show 
clear  picture  for  polarization  of  these  events,  though  MPD  HiV  method  gave  some¬ 
what  lower  standard  deviation  resulting  in  slight  preferrence  for  SV-type  of  waves. 
One  of  the  reasons  for  poor  performance  of  polarization  method  in  these  conditions 
could  be  the  very  rough  edge  of  ice  lead  in  the  vicinity  of  geophone  arrays.  The  failure 


117 


Figure  4-10:  Comparison  of  model  dispersion  and  the  phase  velocity  estimates  ob¬ 
tained  by  conventional  beamforming  of  events  #2  and  #3. 

of  polarization  processing  for  these  events  could  also  help  to  explain  why  the  better 
agreement  between  theory  and  model  for  edge  waves  has  not  been  achieved. 


4.6  Conclusions 

The  described  experiment  [6]  provided  the  evidence  for  the  existence  of  edge  waves, 
propagating  along  a  newly  opened  lead.  The  waves  exhibit  a  quasi-periodic  behavior 
suggesting  some  kind  of  stick-slip  generation  mechanism  somewhere  along  the  length 
of  the  lead.  The  propagation  characteristics  of  these  waves  were  determined  using 
seismic  wavenumber  estimation  techniques.  In  the  low-frequency  limit  the  dispersion 
can  be  modeled  approximately  by  an  interaction  at  the  lead  edges  of  the  lowest 
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order,  antisymmetric  modes  of  the  infinite  plate.  However,  more  modeling  is  required 
to  understand  the  nature  of  this  new  ‘edge  wave’  phenomenon,  in  particular  in  regard 
to  achieving  a  better  match  of  dispersion  behavior. 
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Chapter  5 


Event  Identification 


5.1  Overview 

An  important  component  of  ice  mechanics  processes,  one  which  is  critical  to  our 
understanding,  is  the  type  of  fault  motion  during  an  ice  event.  The  type  of  fault 
motion  is  defined  by  the  relation  between  its  direction  of  fault  motion  and  some 
fixed  coordinate  directions.  Each  type  of  fault  has  its  own  unique  radiation  pattern 
which  can  be  used  to  identify  it.  With  multi-component  data  the  differences  between 
radiation  patterns  for  several  orthogonal  components  of  motion  also  can  be  exploited. 

This  chapter  starts  by  defining  the  fault  orientation  parameters  and  the  three 
modes  of  ice  cracking.  The  simulated  radiation  patterns  for  the  three  orthogonal 
components  of  motion  in  different  types  of  ice  cracks  are  then  discussed,  and  simplified 
analytical  models  for  ice  crack  radiation  patterns  are  introduced.  These  analytical 
models  are  used  in  a  procedure  to  identify  events  which  is  discussed  in  the  next 
section.  The  procedure  includes  the  determination  of  the  directional  characteristics 
of  an  ice  event  from  the  geophone  array  data. 

The  analytical  models  are  fitted  in  turn  to  the  event’s  directional  character¬ 
istics;  the  best-fit  model  gives  the  most  probable  event  mechanism  (or  type  of  fault 
motion  in  the  event).  In  the  next  section  the  results  of  the  determination  of  the  event 
mechanisms  for  several  groups  of  events  are  discussed. 
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North 


Figure  5-1:  Definition  of  the  fault-orientation  parameters  (strike  dip  S)  and  slip- 
direction  (rake  A)  according  to  [1]. 

5.2  Fault  types  definitions 

5.2.1  Fault  orientation  description 

The  definitions  of  the  fault-orientation  parameters  and  slip,  strike  and  dip  directions 
are  given  according  to  [l](p.  106)  on  Fig.  5-1.  On  this  figure  the  strike  angle  is 
measured  clockwise  from  the  north,  with  the  fault  dipping  down  to  the  right  of  the 
strike  direction.  Note  that  0  <  0*  <  The  strike  direction  could  also  be  defined  as 
the  direction  along  the  intersection  of  the  fault  plane  with  the  horizontal  plane  with 


Figure  5-2:  Fault  motion  with  a  dip  angle  S  for  3  different  modes  of  cracking  (ac¬ 
cording  to  [22]):  (a)  tensile  crack,  (b)  dip-slip,  (c)  strike  slip.  On  this  plot  the  x-axis 
corresponds  to  direction  to  the  North,  2:-axis  -  to  the  downward  direction. 
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Figure  5-3:  Sound  profile  in  the  water  on  April  22, 1994.  Solid  line  shows  the  measured 
profile,  the  dashed  one  -  the  approximated  one.  The  difference  between  two  does  not 
exceed  0.5  m/s  which,  I  believe,  is  less  than  time  variability  of  the  profile.  Upper  two 
meters  of  depth  are  occupied  by  the  ice  layer  (compressional  speed  -  3500  m /sec  and 
shear  wave  speed  -  1800  m/sec). 

its  sign  chosen  according  to  the  above  convention.  The  dip  angle  5  is  measured  down 
from  the  horizontal:  0  <  5  <  27r.  It  represents  the  angle  between  the  fault  plane  and 
horizontal  plane. 

5.2.2  Three  modes  of  cracking 

The  fault  motion  for  3  modes  of  cracking  according  to  [22]  is  shown  on  Fig.  5-2. 
On  this  plot  the  a;-axis  points  North,  and  the  2-axis  points  down.  Comparing  with 
FiS-  can  say  that  strike-slip  corresponds  to  a  fault  motion  with  rake  angle 

0  or  180  ,  and  dip-slip  to  A  =  90°.  The  displacements  are  in  the  plane  of  crack 
and  normal  to  the  crack  edge  for  dip-slip,  and  they  are  in  the  fault  plane,  parallel  to 


^  left-lateral  and  right-lateral  strike-slip  faults  respectively 
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the  leading  edge  of  the  crack  in  the  strike-slip  mode.  These  facts  clarify  the  origin 
of  the  names  “strike-slip”  and  “dip-slip”  (one  is  in  the  strike  direction,  the  other 
in  the  dip  direction,  according  to  Fig.  5-1).  Note  that  in  the  seismology  literature 
(see,  for  example,  [1]  p.  106)  the  above  terms  are  reserved  for  cracks  satisfying  the 
additional  condition  that  the  dip  angle  is  5  =  90°.  For  the  tensile  crack  the  rake  angle 
cannot  be  defined,  because  the  particle  displacement  is  normal  to  the  fault  surface 
(see  Fig.  5-2a). 

Sometimes  tensile  mode  is  also  called  Mode  I  or  opening  mode,  dip-slip  - 
Mode  II  (in-plane  shear  mode)  or  sliding  mode,  and  strike-slip  -  Mode  III  (anti-plane 
shear)  or  tearing  mode. 

5.3  Radiation  patterns  for  different  types  of  ice 
cracks 

The  Arctic  environment  used  for  the  calculation  of  radiation  patterns  for  different 
types  of  events  consisted  of  10  layers.  The  first  layer  was  a  vacuum  half-space,  the 
second  one  was  a  homogeneous  elastic  ice  of  2  m  thickness  (according  to  experimental 
records).  The  next  7  layers  were  water  layers  with  a  sound  speed  gradient  for  which 
analytical  solutions  could  be  obtained  in  the  form  of  Airy  functions.  These  water 
layers  were  obtained  as  an  approximation  to  the  measured  sound  speed  profile  on 
April  22,  1994  (see  Fig.  5-3).  The  last  layer  was  a  fluid  halfspace  with  constant 
sound  speed  because  no  information  was  available  about  the  sea  bottom  and  the 
contributions  from  the  bottom  on  the  energy  in  the  ice  could  be  neglected  for  all 
practical  purposes. 

The  calculations  of  the  radiation  patterns  for  different  types  of  ice  events  in  the 
ice  layers  were  made  using  the  0 ASP  3D  module  of  the  seismo-acoustic  computational 
package  OASES  created  by  H.  Schmidt  [44].  For  each  event  type  the  depth  of  the 
source  and  the  receiver  was  assumed  to  be  0.2  m. 
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Figure  5-4:  Radiation  pattern  for  a  dip-slip  crack  with  a  dip  angle  of  5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  vertical  component  of  velocity. 

5.3.1  Shear  crack  radiation  patterns 

The  horizontal  radiation  patterns  for  a  dip-slip  crack  with  a  dip  angle  of  5  =  0°  are 
shown  on  Fig.  5-4  (vertical  component  of  velocity),  Fig.  5-5  (radial  component  of 
velocity),  and  Fig.  5-6  (transverse  component  of  velocity).  The  radiation  patterns  for 
a  dip-slip  crack  with  a  dip  angle  of  =  90°  are  not  shown  because  they  are  practically 
identical  to  those  for  a  dip-slip  crack  with  a  dip  angle  of  5  =  0°.  Three  components  of 
velocity  in  the  near-field  of  the  strike-slip  crack  with  a  dip-angle  of  5  =  0°  are  shown 
on  Figs.  5-7,  5-8,  and  5-9,  respectively.  The  radiation  patterns  for  a  strike-slip  crack 
with  a  dip-angle  of  S  =  90°  are  shown  on  Figs.  5-10,  5-11,  and  5-12,  respectively. 

Note  that  for  a  dip-slip  crack  vertical  and  radial  components  of  velocity  have 
qualitatively  identical  radiation  patterns,  while  the  transverse  component  has  a  radi¬ 
ation  pattern  which  could  be  obtained  from  the  vertical  radiation  pattern  by  rotating 
it  by  90  .  The  same  observation  also  applies  to  a  strike-slip  crack  with  a  dip-angle 
of  5  =  0  .  For  a  strike-slip  crack  with  a  dip-angle  of  5  =  90°,  the  respective  rotation 
angle  should  be  45°. 
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Range  (km) 


Figure  5-5:  Radiation  pattern  for  a  dip-slip  crack  with  a  dip  angle  of  5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  radial  component  of  velocity. 


Range  (km) 

Figure  5-6:  Radiation  pattern  for  a  dip-slip  crack  with  a  dip  angle  of  (5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  transverse  component  of  velocity. 


Figure  5-7:  Radiation  pattern  for  a  strike-slip  crack  with  a  dip  angle  of  5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  vertical  component  of  velocity. 


Figure  5-8:  Radiation  pattern  for  a  strike-slip  crack  with  a  dip  angle  of  5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  radial  component  of  velocity. 
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Figure  5-9:  Radiation  pattern  for  a  strike-slip  crack  with  a  dip  angle  of  5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  transverse  component  of  velocity. 


Range  (km) 


Figure  5-10:  Radiation  pattern  for  a  strike-slip  crack  with  a  dip  angle  of  5  =  90°  i: 
the  near-field  of  the  source  at  a  frequency  of  20  Hz:  vertical  component  of  velocity. 


Range  (km) 

Figure  5-11:  Radiation  pattern  for  a  strike-slip  crack  with  a  dip  angle  of  6  =  90°  in 
the  near-field  of  the  source  at  a  frequency  of  20  Hz:  radial  component  of  velocity. 
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Figure  5-12:  Radiation  pattern  for  a  strike-slip  crack  with  a  dip  angle  of  5  =  90°  in 
the  near-field  of  the  source  at  a  frequency  of  20  Hz:  transverse  component  of  velocity. 
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Figure  5-13:  Radiation  pattern  for  a  tensile  crack  with  a  dip  angle  of  5  =  0°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  vertical  component  of  velocity. 

5.3.2  Tensile  crack  radiation  patterns 

For  a  tensile  crack,  experimentation  with  the  values  of  the  moment  tensor  compo¬ 
nents  showed  that  the  radiation  pattern  is  insensitive  to  the  relative  values  of  these 
components  for  a  fixed  value  of  dip  angle  5.  Calculations  of  radiation  patterns  for 
a  tensile  crack  (for  both  values  of  dip-angle)  were  therefore  carried  out  only  for  the 
case  when  the  only  non-zero  component  of  the  moment  tensor  is  in  the  X-direction. 
For  a  dip-angle  of  5  =  0°  the  radiation  of  a  tensile  crack  is  omnidirectional,  so  only 
the  vertical  component  of  the  ice  velocity  is  shown  (see  Fig.  5-13).  For  a  dip-angle  of 
S  =  90°,  all  three  components  of  velocity  are  shown  on  Figs.  5-14,  5-15,  and  5-16. 

Note  that  while  the  vertical  and  radial  components  of  velocity  for  a  5  =  90° 
tensile  crack  have  qualitatively  identical  radiation  patterns,  the  transverse  component 
has  a  quite  distinct  quadrupole  radiation  pattern  resembling  the  vertical  component 
of  a  strike-slip  crack  with  a  dip-angle  of  5  =  90°. 
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Range  (km) 

Figure  5-14:  Radiation  pattern  for  a  tensile  crack  with  a  dip  angle  of  <5  =  90°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  vertical  component  of  velocity. 
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Figure  5-15:  Radiation  pattern  for  a  tensile  crack  with  a  dip  angle  of  5  =  90°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  radial  component  of  velocity. 
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Figure  5-16:  Radiation  pattern  for  a  tensile  crack  with  a  dip  angle  of  <5  =  90°  in  the 
near-field  of  the  source  at  a  frequency  of  20  Hz:  transverse  component  of  velocity. 

5.3.3  Analytical  models  for  crack  radiation  patterns 

The  radiation  patterns  for  different  modes  of  cracking  can  be  represented  analytically 
by  using  three  fundamental  radiation  patterns  (monopole,  dipole  and  quadrupole) 
and  one  radiation  pattern  obtained  by  the  superposition  of  a  monopole  and  a  dipole 
(further  called  MD-type  radiation  pattern).  As  a  result,  a  dip-slip  crack  with  any  dip- 
angle  and  a  strike-slip  crack  with  a  dip-angle  of  5  =  0°  can  be  represented  analytically 
as  a  dipole;  a  strike-slip  crack  with  a  dip-angle  of  5  =  90°  corresponds  to  a  quadrupole; 
a  tensile  crack  with  a  dip-angle  of  (5  =  0°  radiates  as  a  monopole  source,  while  a  tensile 
crack  with  a  dip-angle  of  <5  =  90°  corresponds  to  a  complex,  MD-type  analytical  model 
of  radiation  pattern.  Analytical  expressions  for  these  four  models  of  radiation  pattern 
are  as  follows: 


const 

for  monopole. 

sin^ 

for  dipole. 

sin  2(p 

for  quadrupole, 

l-\cos2ip 

for  MD-type. 
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Crack  type 

S 

Pattern  type 

Transverse  component 

Dip-slip 

0° 

Dipole 

Rotated  90° 

90° 

Dipole 

Rotated  90° 

Strike-slip 

0° 

Dipole 

Rotated  90° 

90° 

Quadrupole 

Rotated  45° 

Tensile 

0° 

Monopole 

No  rotation 

90° 

MD-type 

No  rotation 

Table  5.1:  Correspondence  between  analytical  models  and  crack  types.  The  last 
column  shows  how  the  radiation  pattern  for  the  transverse  component  can  be  obtained 
from  that  of  the  radial  or  vertical  component. 


We  summarize  the  correspondence  between  analytical  models  and  crack  types 
in  the  Table  5.1. 

5.4  Procedure  for  fitting  analytical  radiation  pat¬ 
terns  to  experimental  data 

5.4.1  Model  parameterization 

The  comparison  between  experimental  data  (relative  signal  levels  on  different  geo¬ 
phones)  and  analytical  radiation  patterns  has  been  made  separately  for  each  of  the 
components  of  motion.  The  comparison  has  been  carried  out  for  levels  in  dB  both 
for  experimental  data  and  radiation  patterns.  There  is  a  slight  difficulty  here  be¬ 
cause,  as  one  can  see  from  Equation  (5.1),  some  of  the  analytic  patterns  (dipole  and 
quadrupole)  have  zeros  which  give  rise  to  infinite  levels  in  dB.  To  avoid  this  problem 
one  needs  to  specify  the  dynamic  range  of  the  radiation  patterns  beforehand.  It  is 
convenient  to  specify  such  a  dynamic  range  from  0  dB  to,  say,  L^ax  dB.  In  L^ax  we 
have  a  parameter  of  radiation  pattern  that  could  be  changed  to  adjust  the  quality 
of  fit  between  the  analytical  radiation  pattern  and  experimental  data.  The  other 
such  parameter  is  the  angle  of  rotation  of  the  original  radiation  pattern  given  by 
Equation  (5.1).  It  is  obvious  that  this  parameter  is  not  useful  for  monopole-style 
radiation  patterns  due  to  the  omnidirectional  character  of  that  pattern.  For  other 
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radiation  patterns  to  achieve  rotation  by  (po  radians  in  the  clockwise  direction  we 
need  to  substitute  (p  (p  +  (po  in  Equation  (5.1). 

As  a  result,  we  have  two  parameters  for  adjusting  the  fit  between  the  ana¬ 
lytical  radiation  pattern  and  the  experimentally  determined  signal  levels  on  different 
geophones: 

Lmax  -  maximum  level  of  radiation  pattern  in  dB; 

(po  -  angle  of  rotation  of  original  analytical  radiation  pattern  given  by  Equation  (5.1). 
According  to  Fig.  5-1,  (po  corresponds  to  the  strike  angle. 

5.4.2  Preprocessing  of  experimental  data  for  determining  the 
radiation  patterns 

Now,  let  us  consider  experimental  data.  First,  the  event  is  located  in  space.  Knowing 
its  location,  the  respective  direction  to  each  geophone  is  easily  determined.  Denoting 
Cartesian  coordinates  of  the  event  as  {xo,yo)  and  the  ith  geophone  -  as  {xi,yi), 
we  have  the  polar  coordinates  of  the  ith.  geophone  in  the  source  coordinate  system 
(centered  at  the  event  location): 

Ri  =  (5.2) 

ipi  =  arctan  — — (5-3) 

Xi  -  Xo 

The  next  step  is  to  determine  the  radiation  level  in  the  direction  given  by  (pi  from 
Equation  (5.3).  This  level  for  each  component  of  motion  is  determined  as  the  peak 
intensity  level  in  the  considered  arrival:  Lj.  As  a  result,  for  for  each  component  of 
motion  we  have  5  points  ^  given  by  their  polar  coordinates  (Lj,  cpi)  which  determine 
(sometimes  ambiguously)  the  radiation  pattern  of  the  event  in  question. 


^  Each  RLAM  unit  had  only  5  geophones 
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5.4.3  Model  pattern  fitting  to  experimental  points 

Now,  again  for  each  component  of  motion  separately,  all  available  analytical  models  of 
radiation  patterns  (see  Equation  (5.1))  are  fitted  to  the  experimental  data  as  specified 
above.  For  each  particular  model  the  best-fit  and  are  found.  The  qualities 
of  fit  for  different  models  are  compared  and  the  best-fit  model  Aio  is  chosen,  having 
corresponding  best-fit  parameters  and  The  quality  of  fit  for  model  Mj 

is  determined  by  the  standard  deviation  aj>^.  of  the  experimental  points  from  the 
model  radiation  pattern.  The  deviation  of  individual  experimental  points  from  the 
analytical  model  is  determined  as  a  minimal  distance  from  the  experimental  points 
to  the  curve  representing  the  analytical  model  in  polar  coordinates.  Due  to  the 
non-linear  closed-form  expressions  for  most  analytical  models  (see  Equation  (5.1)) 
the  least-square  fit  for  individual  experimental  points  (in  order  to  find  the  minimal 
distance  from  the  point  to  analytical  curve)  was  carried  out  numerically  by  sampling 
the  analytical  curve  in  polar  coordinates  with  sufficient  angular  resolution  and  finding 
the  point  in  this  discretized  set  with  minimal  distance  to  the  experimental  point  in 
question. 

In  other  words,  we  have  a  set  of  models  £  M  each  parameterized  by  T^ax 
and  (fo-  The  best-fit  model  can  be  expressed  as 

A1„  =  argmto  (5.4) 

The  flow-chart  for  the  model  fit  for  each  component  of  motion  is  shown  on  Fig.  5-17. 
Essentially,  this  procedure  is  a  simple  example  of  functional  (or  variational)  analysis. 
The  last  stage  of  the  fit  process  is  to  compare  best-fit  models  for  different  compo¬ 
nents  of  motion.  According  to  Table  5.1,  we  should  get  identical  results  (in  terms  of 
model  and  its  best-fit  parameters)  for  radial  and  vertical  components,  whereas  for  the 
transverse  component,  tpo  should  differ  by  45°  for  a  strike-slip  fault  with  a  dip  angle 
of  90  and  by  90°  for  a  strike-slip  fault  with  a  dip  angle  of  0°  and  for  a  dip-slip  fault 
with  a  dip  angle  of  0°  or  90°.  If  we  observe  the  mentioned  above  differences  in  best-fit 
parameters  for  different  components  of  motion  (for,  at  least,  two  components),  this 
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Model  input 


Figure  5-17:  Flow-chart  of  the  analytical  radiation  pattern  fit  for  each  component  of 
motion. 

should  increase  our  confidence  in  the  determination  of  the  source  mechanism. 

5.5  Qualitative  determination  of  event  mechanisms 

Using  the  above  correspondence  between  crack  types  and  the  radiation  pattern  mod¬ 
els,  a  qualitative  determination  of  the  source  mechanisms  of  ice  events  was  made  for 
several  groups  of  events.  In  the  results  discussed  below  the  standard  deviation  of  the 
best-fit  model  from  the  experimental  points  is  presented  on  a  linear  scale  as  a  fraction 
of  the  amplitude  of  the  mainlobe  of  the  model  (expressed  in  %). 
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Figure  5-18:  Locations  of  five  events  in  Cluster  #1.  Triangle  markers  are  used  to 
show  the  locations  of  geophones,  while  filled  circles  represent  the  locations  of  events. 
The  numbers  near  the  circle  markers  show  approximate  arrival  times  in  seconds  since 
the  beginning  of  tape  RLAM-26  for  corresponding  events. 

5.5.1  Cluster  #1 

The  first  group  was  Cluster  #1,  a  group  of  five  events  (for  the  locations  of  the  events 
see  Fig.  5-18)  spanning  approximately  10  sec  at  the  beginning  of  tape  RLAM-26  with 
very  similar  time-domain  behavior.  The  time  series  of  the  horizontal  component  of 
velocity  for  four  events  from  that  cluster  are  shown  on  Fig.  5-19.  One  can  observe 
that  for  all  four  events  three  geophones  (2-4)  have  a  similar  arrival  pattern,  while 
geophones  #1  and  #5  have  quite  different  arrival  patterns.  The  locations  of  the 
events  for  this  cluster  were  determined  by  a  standard  least-square  fit  to  the  arrival 
times. 

For  event  #5  (approximately  39  sec  from  the  beginning  of  tape  RLAM-26) 
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Figure  5-20:  Best-fit  radiation  pattern  models  for  event  #5  from  Cluster  #1  (ap¬ 
proximately  39  sec  from  the  beginning  of  tape  RLAM-26):  (a)  vertical  component 
of  motion,  cpo  =  —39°,  =  4.5%;  (b)  radial  component  of  motion  (to  source), 

^0  =  —35°,  gm  =  3.1%.  Red  circle  markers  correspond  to  experimental  data  from 
the  five  geophones. 


Figure  5-21:  Best-fit  radiation  pattern  models  for  event  #7  from  Cluster  #1  (ap¬ 
proximately  46  sec  from  the  beginning  of  tape  RLAM-26):  (a)  vertical  component  of 
motion,  =  41°,  gm  =  4.5%;  (b)  radial  component  of  motion  (to  source),  Pq  —  45°, 
Cm  =  2.4%.  Red  circle  markers  correspond  to  experimental  data  from  the  five  geo¬ 
phones. 
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the  above  procedure  yielded  a  quadrupole  as  the  best-fit  model,  which  corresponds  to 
a  strike-slip  fault  with  a  dip  angle  of  90°.  In  Fig.  5-20  the  best-fit  models  are  shown 
for  the  vertical  and  radial  components.  These  models  gave  quite  close  values  of  strike 
angles  (-39°  and  -35°,  respectively).  The  fit  for  the  transverse  component  in  this 
case  is  not  acceptable  and  therefore  is  not  shown  here.  Similar  results  for  event  #7 
(approximately  46  sec  from  the  beginning  of  tape  RLAM-26)  are  shown  on  Fig.  5-21. 
Again,  the  strike-slip  fault  mechanism  (with  a  dip  angle  of  90°)  is  suggested  for  the 
event.  The  values  of  strike  angles  in  this  case  are  41°  for  the  vertical  component,  and 
45°  for  the  radial  component. 
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Figure  5-22:  Locations  of  seven  events  at  the  beginning  of  tape  RLAM-26.  Triangle 
markers  are  used  to  show  the  locations  of  geophones,  while  filled  circles  represent  the 
locations  of  events.  The  numbers  near  circle  markers  show  approximate  arrival  times 
in  seconds  since  the  beginning  of  tape  RLAM-26  for  corresponding  events. 

5.5.2  Other  events  at  the  beginning  of  tape  RLAM-26 

Several  other  events  at  the  beginning  of  tape  RLAM-26  were  also  successfully  located 
by  standard  least-square  fit  to  arrival  times.  Their  locations  are  shown  on  Fig.  5-22. 

For  these  events,  an  especially  good  fit  was  obtained  for  event  #10  (approx¬ 
imately  64  sec  from  the  beginning  of  tape  RLAM-26).  The  discrimination  between 
different  models  was  also  outstanding.  The  best  fit  was  achieved  for  a  quadrupole  radi¬ 
ation  pattern,  corresponding  to  a  strike-slip  fault  with  a  dip  angle  of  90°.  On  Fig.  5-23 
the  best-fit  models  are  shown  for  the  vertical  and  radial  components.  The  quality 
of  fit  is  better  appreciated  by  comparing  with  fits  by  two  other  models  (dipole  and 


140 


90 


15 


90 


20 


270  270 

a)  b) 

Figure  5-23:  Best-fit  models  of  radiation  patterns  for  event  #10  (approximately  64  sec 
from  the  beginning  of  tape  RLAM-26):  (a)  vertical  component  of  motion,  (pQ  =  16°, 
om  —  4.4%;  (b)  radial  component  of  motion  (to  source),  Pq  =  —1°,  gm  =  2.1%.  Red 
circle  markers  correspond  to  experimental  data  from  the  five  geophones. 
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a)  b) 

Figure  5-24:  Results  of  fits  by  two  models  (dipole  and  MD-type)  for  the  vertical 
component  of  motion  during  event  #10  (approximately  64  sec  from  the  beginning 
of  tape  RLAM-26):  (a)  dipole  model,  ipQ  =  60°,  gm  =  10.6%;  (b)  MD-type  model, 
ipQ  =  60°,  gm  =  15.0%.  Red  circle  markers  correspond  to  experimental  data  from  the 
five  geophones. 
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Figure  5-25:  Results  of  fits  by  two  models  (monopole  and  dipole)  for  the  vertical 
component  of  motion  during  event  #16  (approximately  168  sec  from  the  beginning 
of  tape  RLAM-26):  (a)  monopole  model,  om  =  19.2%;  (b)  dipole  model,  = 
—20°,  gm  =  4.1%.  Red  circle  markers  correspond  to  experimental  data  from  the  five 
geophones. 
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a)  b) 

Figure  5-26:  Results  of  fits  by  two  models  (quadrupole  and  MD-type)  for  vertical 
component  of  motion  during  event  #16  (approximately  168  sec  from  the  beginning 
of  tape  RLAM-26):  (a)  quadrupole  model,  =  —10°,  gm  =  4.1%;  (b)  MD-type 
model,  ^0  =  —40°,  gm  =  11-0%.  Red  circle  markers  correspond  to  experimental  data 
from  the  five  geophones. 
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MD-type).  The  results  of  fits  for  these  models  are  shown  on  Fig.  5-24.  A  monopole 
gives  an  even  worse  performance  than  an  MD-type  radiation  pattern  (the  deviation 
from  the  monopole  is  two  times  larger  than  that  from  the  MD-type  pattern).  Such 
a  good  discrimination  between  different  models  was  possible  in  this  case  because  the 
event  was  located  in  the  immediate  vicinity  of  the  receiving  array  (see  Fig.  5-22). 

Event  ^^16  (approximately  168  sec  from  the  beginning  of  tape  RLAM-26)  rep¬ 
resents  just  the  opposite  case.  It  is  located  quite  far  (about  400  m)  from  the  geophone 
array  (see  Fig.  5-22),  and  the  experimental  data  have  a  small  angular  spread.  As  a 
result,  we  could  not  really  discriminate  between  the  different  models  and  therefore 
could  not  assign  a  corresponding  event  mechanism  (compare  Fig.  5-25  and  Fig.  5-26). 
This  case  shows  one  of  the  reasons  why  sometimes  even  a  qualitative  determination 
of  the  event  mechanism  cannot  be  made  for  a  successfully  located  event. 

5.5.3  Events  at  the  beginning  of  file  #5  of  tape  RLAM-26 

The  last  group  of  events  for  which  the  event  mechanism  determination  has  been 
carried  out  was  the  group  of  events  at  the  beginning  of  file  #5  on  tape  RLAM-26. 
These  events  were  located  using  the  MPD  method  described  elsewhere  in  this  thesis. 
The  locations  of  the  events  discussed  below  are  shown  on  Fig.  5-27. 

The  best-fit  radiation  pattern  models  for  event  #5  (approximately  165  sec 
from  the  beginning  of  file  #5  on  tape  RLAM-26)  are  shown  Fig.  5-28  (vertical  com¬ 
ponent  of  motion)  and  Fig.  5-29  (radial  and  transverse  components  of  motion).  For 
all  components  the  best  fit  is  obtained  by  a  quadrupole  model  which  implies  a  strike- 
slip  fault  with  a  dip  angle  of  90°.  This  conclusion  is  also  supported  by  the  value  of 
the  difference  in  the  best-fit  parameter  between  vertical  and  transverse  compo¬ 
nents,  38°,  and  between  radial  and  transverse  components  -  39°.  In  both  cases  this 
difference  is  quite  close  to  45°  as  it  should  be  according  to  numerical  simulations  (see 
Table  5.1). 

The  difference  between  (pQ  for  the  radial  and  transverse  components  is  even 
closer  to  45°  for  event  #11  (approximately  374  sec  from  the  beginning  of  file  #5 
on  tape  RLAM-26).  For  both  components  the  best  fit  is  achieved  by  a  quadrupole 
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Figure  5-27:  Locations  of  four  events  from  file  ^5  of  tape  RLAM-26.  Triangle  markers 
are  used  to  show  the  locations  of  geophones,  while  filled  circles  show  the  locations 
of  events.  The  numbers  near  the  circle  markers  show  approximate  arrival  times  in 
seconds  since  the  beginning  of  file  #5  of  tape  RLAM-26  for  corresponding  events. 
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Figure  5-28:  Best-fit  model  of  radiation  pattern  for  the  vertical  component  of  motion 
during  event  #5  (approximately  165  sec  from  the  beginning  of  file  #5  on  tape  RLAM- 
26).  (fo  =  8°,  aM  =  8.2%.  Red  circle  markers  correspond  to  experimental  data  from 
the  five  geophones. 


Figure  5-29:  Best-fit  models  of  radiation  patterns  for  event  #5  (approximately  165  sec 
from  the  beginning  of  file  #5  on  tape  RL AM-26):  (a)  radial  component  of  motion  (to 
source),  (po  =  9°,  aM  =  10.2%;  (b)  transverse  component  of  motion,  (fo  =  -30°,  aM  = 
9.8%.  Red  circle  markers  correspond  to  experimental  data  from  the  five  geophones. 
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a)  b) 


Figure  5-30:  Best-fit  models  of  radiation  patterns  for  event  #11  (approximately 
374  sec  from  the  beginning  of  file  #5  on  tape  RLAM-26):  a)  radial  component 
of  motion  (to  source),  ipQ  =  6°,  gm  =  6.9%;  b)  transverse  component  of  motion, 
^0  =  —38°,  gm  =  6.3%.  Red  circle  markers  correspond  to  experimental  data  from 
the  five  geophones. 


Figure  5-31:  Best-fit  models  of  radiation  patterns  for  event  #12  (approximately 
376  sec  from  the  beginning  of  file  #5  on  tape  RLAM-26):  a)  vertical  component 
of  motion,  =  —31°,  gm  =  5.5%;  b)  transverse  component  of  motion,  =  14°, 
—  7.2%.  Red  circle  markers  correspond  to  experimental  data  from  the  five  geo¬ 
phones. 
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Figure  5-32:  Best-fit  models  of  radiation  patterns  for  event  #9  (approximately  361  sec 
from  the  beginning  of  file  #5  on  tape  RL AM-26):  a)  radial  component  of  motion 
(to  source),  (pQ  =  -34°,  gm  =  10.9%;  (b)  transverse  component  of  motion,  ipQ  = 
27°,  gm  =  2.3%.  Red  circle  markers  correspond  to  experimental  data  from  the  five 
geophones. 


radiation  pattern  (see  Fig.  5-30)  which  corresponds  to  a  strike-slip  fault  with  a  dip 
angle  of  90°.  The  same  event  mechanism  determination  was  made  for  event  ^12  (see 
Fig.  5-31).  In  that  case  the  difference  in  pq  was  exactly  45°  which  gives  additional 
assurance  as  to  a  correct  event  mechanism  identification. 

For  event  #9  (approximately  361  sec  from  the  beginning  of  file  #5  on  tape 
RLAM-26)  the  difference  in  pq  between  the  radial  and  transverse  components,  61° 
(see  Fig.  5-32),  was  not  as  close  to  45°  as  for  events  #5,  #11  and  #12,  but  still 
within  reasonable  bounds  (61°  is  almost  two  times  closer  to  45°  than  to  90°).  The 
event  mechanism  in  this  case  is  identified  as  a  strike-slip  fault  with  a  dip  angle  of  90°. 
For  this  event  the  standard  deviation  of  the  quadrupole  model  of  radiation  pattern 
as  a  function  of  two  parameters,  strike  angle  and  -Tmax,  is  shown  on  Fig.  5-33. 
Different  colors  of  contours  on  this  figure  correspond  to  different  levels  of  g  with  blue 
being  the  lowest  level  and  red  the  highest.  The  minimum  (or  best  fit)  is  achieved  for 
Pq  =  27°  and  I<mav=14  dB. 

The  results  of  the  determination  of  fault  type  for  several  ice  events  are  shown 
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Figure  5-33:  Standard  deviation  between  experimental  data  and  the  quadrupole  nodes 
as  a  function  of  strike  angle  yjo  and  -Lmax  (designated  as  “Min.  level”  on  plot),  for 
transverse  component  of  motion  during  event  #9  (approximately  361  sec  from  the 
beginning  of  file  on  tape  RL AM-26).  Different  colors  of  contours  correspond  to 
different  levels  of  a  with  blue  being  the  lowest  level  and  red  the  highest.  The  minimum 
(or  best  fit)  is  achieved  for  (po  =  27°  and  Tmax=14  dB. 
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Figure  5-34:  Results  of  the  determination  of  fault  type  for  several  ice  events  from  “ice 
island”  site. 
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Figure  5-35:  Directions  of  the  fault  motion  for  several  events  from  “ice  island”  site. 


on  Fig.  5-34.  These  results  indicate  that  strike-slip  faults  with  dip  angles  of  90°  are 
the  most  probable  source  mechanism.  Fig.  5-35  depicts  direction  of  the  fault  motion 
for  ice  events  discussed  in  this  Chapter. 


5.6  Conclusions 

In  this  chapter  ice  event  identification  was  successfully  carried  out  for  several  groups 
of  events  by  comparing  the  quality  of  fit  of  different  analytical  radiation  patterns 
to  experimental  data.  The  results  of  the  event  identification  suggest  that  strike-slip 
faults  with  dip  angles  of  90°  are  the  most  probable  source  mechanism.  Additional 
evidence  in  favor  of  this  conclusion  is  supplied  by  the  differences  in  orientation  of  the 
best-fit  radiation  patterns  for  orthogonal  components  of  motion  which  are  very  close 
to  the  values  given  by  numerical  models. 

This  choice  of  strike-slip  fault  with  a  dip  angle  of  90°  as  the  most  probable 
source  mechanism  for  processed  ice  events  is  somewhat  to  the  contrary  of  J.S.  Kim’s 
conclusion  [22]  which  was  based  on  the  relative  smallness  of  the  radiation  from  longi¬ 
tudinal  waves  compared  to  the  acoustic  mode.  He  suggested  that  dip-slip  craeks  with 
dip  angles  of  5  «  0°  and  90°  and  strike-slip  cracks  with  a  dip  angle  of  5  «  0°  are  the 
most  probable  source  mechanisms.  I  think  this  discrepancy  may  arise  from  the  fact 
that  most  of  the  events  processed  here  do  not  radiate  efficiently  into  the  water. 
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Chapter  6 


Event  parameter  estimation 

6.1  Overview 

The  similarity  between  physical  processes  in  earthquakes  and  ice  cracking  has  already 
been  noted  in  the  introduction.  Based  on  that  similarity  we  will  use  established 
methods  from  earthquake  seismology  to  estimate  the  two  most  important  parameters 
of  faulting  processes:  stress  drop  and  seismic  moment.  These  methods  can  be  divided 
into  two  broad  classes:  spectral  and  slip  history.  The  spectral  method  of  parameters 
estimation  will  be  discussed  in  the  first  section  along  with  some  examples  of  data 
processing.  In  the  next  section  the  determination  of  slip  history  will  be  described  and 
the  data  processing  results  (including  stick-slip  motion)  will  be  discussed.  The  chapter 
will  continue  with  the  definition  of  stress  drop  and  seismic  moment  for  ice  events.  Next 
the  procedure  for  the  determination  of  event  parameters  from  the  combined  results  of 
spectral  and  slip  history  methods  will  be  explained.  Finally,  the  results  of  the  event 
parameter  estimation  will  be  discussed  and  a  comparison  with  previous  estimates  will 
be  presented. 
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Figure  6-1:  Standard  model  of  the  displacement  spectrum  of  an  earthquake.  On 
this  figure  S  (the  displacement  spectrum  level  in  dB)  is  plotted  versus  the  logarithm 
of  frequency.  Two  asymptotes  are  usually  associated  with  this  spectrum:  a  high- 
frequency  asymptote  and  a  low-frequency  asymptote.  Their  intersection  defines  the 
corner  frequency  fo- 

6.2  Frequency  domain  processing  for  event  param¬ 
eter  estimation 

6.2.1  Standard  earthquake  model  of  the  displacement  spec¬ 
trum 

In  seismology  the  following  model  of  displacement  spectrum  for  a  plane  fault  is  used 
[48], [28]: 

1.  In  the  high  frequency  limit  the  displacement  spectrum  approaches  an  asymptote 
of  the  form  oo~'^  (usually,  in  seismology  7  is  taken  to  be  2). 

2.  The  intersection  of  this  high-frequency  asymptote  with  the  low  frequency  level 
gives  the  corner  frequency  fo  (see  Fig.  6-1)  which,  according  to  [3],  is  inversely 
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Figure  6-2:  Vertical  displacement  spectra  of  the  event  #3  (approximately  36.5  sec 
from  the  beginning  of  the  tape  RL AM-26).  Different  colors  correspond  to  different 
geophones.  The  high-frequency  roll-off  has  a  trend  (one  of  the  black  lines).  The 
other  black  line  shows  the  low-frequency  trend.  The  intersection  of  the  two  trends 
gives  the  corner  frequency  /o  as  7.8  Hz. 

related  to  the  dimensions  of  the  source  (or  source  duration  T):  /o  =  I/ttT. 

6.2.2  Results  of  corner  frequency  estimation 

The  results  of  corner  frequency  estimation  distinctly  fall  into  two  groups:  the  first 
one  having  only  one  high  frequency  trend  (w“®)  and  the  second  one  having  both 
and  high  frequency  trends. 

The  first  group  is  shown  on  Figs.  6-2,  6-3  and  6-4.  On  these  figures  the  vertical 
displacement  spectra  for  several  ice  events  are  shown  for  all  five  geophones.  On  Fig.  6- 
2  the  results  of  the  processing  of  event  ^3  (approximately  36.5  sec  from  the  beginning 
of  tape  RLAM-26)  are  shown.  For  this  event  the  corner  frequency  was  estimated  to  be 
/o  =  7.8  Hz.  For  event  ^5,  which  occured  approximately  39  sec  from  the  beginning 
of  tape  RLAM-26,  the  spectral  method  gave  the  following  value  of  corner  frequency: 
/o  =  8.8  Hz  (see  Fig.  6-3).  The  results  of  the  processing  of  event  #9  (approximately 
58  sec  from  the  beginning  of  tape  RLAM-26)  are  shown  on  Fig.  6-4.  In  this  case  the 
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Figure  6-3;  Vertical  displacement  spectra  of  event  #5  (approximately  39  sec  from  the 
beginning  of  tape  RLAM-26).  Different  colors  correspond  to  different  geophones.  The 
high-frequency  roll-off  follows  a  (jJ~^  trend  (one  of  the  black  lines).  The  other  black 
line  shows  the  low-frequency  trend.  Their  intersection  gives  the  corner  frequency 
/o  =  8.8  Hz. 


Figure  6-4:  Vertical  displacement  spectra  of  event  #9  (approximately  58  sec  from  the 
beginning  of  tape  RLAM-26).  Different  colors  correspond  to  different  geophones.  The 
high-frequency  roll-off  follows  a  trend  (one  of  the  black  lines).  The  other  black 
line  shows  the  low-frequency  trend.  Their  intersection  gives  the  corner  frequency 
/o  =  8.4  Hz. 
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Figure  6-5:  Vertical  displacement  spectra  of  the  event  #11  (approximately  74  sec  from 
the  beginning  of  tape  RL AM-26).  Different  colors  correspond  to  different  geophones. 
For  this  event  two  high-frequency  trends  can  be  observed.  First,  there  is  a  trend 
(one  of  the  black  lines),  then  a  uj~^  trend  (the  second  black  line).  No  corner  frequency 
could  be  identified  in  this  case. 


Figure  6-6:  Vertical  displacement  spectra  of  the  event  #12  (approximately  80  sec 
from  the  beginning  of  the  tape  RL  AM-26).  Different  colors  correspond  to  different 
geophones.  For  this  event  two  high-frequency  trends  can  be  observed.  First,  there  is 
aa;-3  trend  (one  of  the  black  lines),  then  a  w  ^  trend  (the  second  black  line).  No 
corner  frequency  could  be  identified  in  this  case. 
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Figure  6-7:  Vertical  displacement  spectra  of  event  #13  (approximately  118  sec  from 
the  beginning  of  tape  RL AM-26).  Different  colors  correspond  to  different  geophones. 
For  this  event  two  high-frequency  trends  can  be  observed.  First,  there  is  a  trend 
(one  of  the  black  lines),  then  a  trend  (the  second  black  line).  No  corner  frequency 
could  be  identified  in  this  case. 


corner  frequency  value  was  estimated  to  be  fo  =  8.4  Hz. 

The  second  group  of  events  is  shown  on  Figs.  6-5,  6-6  and  6-7.  The  results 
of  the  processing  for  event  #11  (approximately  74  sec  from  the  beginning  of  tape 
RLAM-26)  are  shown  on  Fig.  6-5,  for  event  #12  (approximately  80  sec  from  the 
beginning  of  the  tape  RLAM-26)  on  Fig.  6-6,  and  for  event  #13  (approximately 
118  sec  from  the  beginning  of  the  tape  RLAM-26)  on  Fig.  6-7.  On  these  figures  the 
vertical  displacement  spectra  for  several  ice  events  are  shown  for  all  five  geophones. 
The  common  feature  of  the  results  of  processing  in  this  group  is  that  the  displacement 
spectra  have  two  high-frequency  trends  starting  with  u}~^  trend  and  ending  with 
trend.  Due  to  the  ambiguity  of  this  situation,  no  corner  frequency  could  be  identified 
in  this  case. 

In  earthquake  mechanics  a  u}~^  trend  is  usually  associated  with  the  process 
of  nucleation  of  the  fault  [1].  On  the  basis  of  only  limited  geophone  coverage  of 
the  faulting  processes  this  hypothesis  could  be  neither  confirmed  nor  denied  for  the 
discussed  ice  events,  but  it  could  be  considered  as  a  plausible  explanation  for  this 
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phenomenon. 


6.3  Slip  history  determination 

6.3.1  Two  models  for  the  slip  history  of  the  fault 

The  two  most  widely  accepted  theories  in  seismology  of  the  slip  history  on  the  fault 
are  models  of  Haskell  [16]  and  Brune  [3]. 

In  the  Haskell’s  model  [16]  (see  Fig.  6-8a)  the  displacement  on  the  fault  u{t) 
could  be  expressed  as 


u{t)  =  Ucx,G{t  —  x/Vr), 


(6.1) 


where  Vr  is  the  rupture  propagation  velocity,  and  G{t)  is  a  ramp  function.  The  ramp 
function  is  zero  at  <  <  0  and  increases  linearly  with  time  until  it  reaches  1  at  t  =  Tr, 
which  is  called  the  rise  time  (for  asymptotic  curves,  the  time  to  reach  90%  of  the 
final  level  of  the  curve  is  sometimes  called  the  rise  time). 

In  Brune’s  model  [3]  (see  Fig.  6-8b)  we  have  for  the  displacement 

u{t)  =  Woo  ^1  —  e“^  ^  for  t  >  0.  (6.2) 

For  the  initial  particle  velocity  this  theory  gives 

u 

Uq  —  Cg 


where: 

a  is  the  initial  stress  level, 

p,  is  the  shear  strength  of  the  medium, 

Cg  (or  as  it  is  usually  denoted  in  seismological  literature)  is  the  shear  wave  speed  of 
the  medium. 
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Figure  6-8:  Two  standard  models  of  slip  history  accepted  in  seismology:  (a)  Haskell’s 
model;  (b)  Brune’s  model.  On  this  figure  for  both  models  the  following  parameters 
are  defined:  uo  —  initial  velocity  (in  text,  it  is  also  called  ito),  Uoo  —  maximum  slip, 
Tr  —  rise  time. 


In  the  seismological  community  Brune’s  model  is  usually  preferred  over  Haskell’s 
model  as  it  more  closely  describes  the  physics  of  faulting.  Brune’s  model  takes  into 
the  account  the  finiteness  of  the  fault  dimensions,  while  Haskell’s  does  not.  For  the 
following  discussion,  Brune’s  model  will  be  used. 

According  to  Kasahara  ([21],  p.  108),  an  average  particle  velocity  associated 
with  faulting  may  be  derived  simply  by  dividing  the  slip  amplitude  by  rise  time.  In 
other  words  we  have 


(6.4) 


where: 

iiQ  is  the  initial  (or  average)  velocity, 

Uoo  is  the  slip  amplitude. 

For  event  parameter  estimation  the  most  important  parameter  is  the  average 
slip  u  which  represents  a  slip  amplitude  Uoo  which  has  been  averaged  in  some  way. 
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Figure  6-9:  Slip  history  of  event  #1  (approximately  27  sec  from  the  beginning  of  tape 
RLAM-26).  Only  the  horizontal  X-component  (true  North)  of  the  displacement  is 
shown  on  this  plot.  Four  thin  lines  correspond  to  4  geophones  appropriately  shifted 
in  time  in  order  to  align  their  peaks.  The  thick  line  shows  the  average  across  4 
geophones  aligned  in  time.  This  curve  gives  a  rise  time  of  approximately  62  msec  and 
a  slip  amplitude  of  5.13  x  10~®  m. 

6.3.2  Results  of  slip  history  determination 

Simple  slip  histories 

On  Fig.  6-9  you  can  see  the  slip  history  of  event  #1  (approximately  27  sec  from 
the  beginning  of  tape  RLAM-26).  Only  the  horizontal  X-component  (true  North)  of 
the  displacement  is  shown  on  this  plot.  Note  that  at  the  end  the  displacement  does 
not  come  to  a  constant  level  as  in  the  model  of  slip  history.  The  oscillations  at  the 
end  of  the  event  could  be  explained  by  hi-pass  filtering  at  the  geophone  itself:  the 
low  frequency  cutoff  of  each  geophone  was  4.5  Hz  at  which  frequency  the  frequency 
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Slip  history  for  2  events 


Figure  6-10:  Slip  history  of  two  events  at  about  38  sec.  from  the  beginning  of  the 
tape  RLAM-26.  All  three  components  of  displacement  are  shown.  On  this  plot  the 
units  of  displacement  are  arbitrary.  The  scales  are  shown  for  separate  events. 
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Figure  6-11:  Slip  history  of  the  first  event  from  Fig.  6-10.  Only  the  horizontal  compo¬ 
nents  of  displacement  are  shown.  On  this  plot  the  units  of  displacement  are  arbitrary. 
The  scale  is  shown  on  the  plot  of  the  Y-component.  The  thick  line  shows  the  average 
across  all  5  geophones  (the  alignment  was  not  really  needed  in  this  case  unlike  the 
case  of  Fig.  6-9). 
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Figure  6-12:  Slip  history  of  the  second  event  from  Fig.  6-10.  Only  the  horizontal 
components  of  displacement  are  shown.  On  this  plot  the  units  of  displacement  are 
arbitrary.  The  scale  is  shown  on  the  plot  of  the  X-component.  The  thick  line  shows 
the  average  across  all  5  geophones  (the  alignment  was  not  really  needed  in  this  case 
unlike  in  Fig.  6-9). 
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response  was  already  3  dB  down.  Looking  at  the  figure  we  could  see  that  the  frequency 
of  the  final  oscillation  is  close  to  the  low  frequency  cutoff. 

The  average  across  the  geophones  in  this  case  gives  a  rise  time  of  62  msec  and 
a  slip  amplitude  of  5.13  x  10“*  m.  From  Equation  (6.4)  we  now  obtain 

Uo  =  8.28  X  10~®  m/s. 

On  Fig.  6-10  the  slip  histories  of  two  events  (approximately  38  sec  from  the 
beginning  of  tape  RLAM-26)  are  shown  for  three  components  of  motion.  For  the 
first  event  (see  Fig.  6-11)  the  determination  of  the  rise  time  could  best  be  achieved 
from  the  Y-component  of  the  displacement  (true  West).  The  rise  time  in  this  case  is 
about  110  msec,  the  slip  amplitude  (before  correction)  is  2.43  x  10“^  m.  According 
to  Equation  (6.4)  that  gives  the  following  value  of  initial  velocity 

uq  =  2.21  X  10“®  m/s. 

For  the  second  event  (see  Fig.  6-12)  the  determination  of  the  rise  time  can  best 
be  achieved  from  the  X-component  of  displacement  (true  North).  The  rise  time  in 
this  case  is  about  120  msec  and  the  slip  amplitude  (before  correction)  is  1.42  x  10“^  m. 
Now,  according  to  Equation  (6.4)  we  get  the  following  value  for  initial  velocity: 

iio  =  1.18  X  10~®  m/s. 

Note  that  above  observation  about  final  oscillations  applies  also  in  these  last 
two  cases. 

Stick-slip  history 

Brace  and  Byerlee  showed  [2]  that  previously  faulted  surfaces  of  rock  may  undergo  a 
series  of  intermittent  small  slips  under  continuing  pressure.  They  called  this  kind  of 
motion  ”  stick-slip”  by  analogy  with  the  term  used  in  materials  science.  Processes  in 
stick-slip  were  described  by  Byerlee  as  follows 


165 


0.03 


>  0.015 

o  0.01 
o 

(D 

>  0.005 

cd 

o 

r  0 

(D 

> 

--0.005 


-0.015 

O.f 


O  0.015 

O  0.01 

o 

a> 

^  0.005 
(d 

o 

> 


Time,  sec. 


Time,  sec. 

(b) 

Figure  6-13:  Time  series  of  the  vertical  component  of  velocity  on  the  geophone  #2 
at  the  beginning  of  tape  RLAM-26.  Note  several  pulses  similar  in  shape  at  approxi¬ 
mately  1.1  sec,  1.4  sec  and  1.8  sec  on  the  top  plot  (a).  The  bottom  plot  (b)  shows  the 
details  of  the  pulse  at  1.4  sec.  The  characteristic  arrival  pattern  indicates  a  flexural 


Figure  6-14:  Autocorrelation  function  for  the  time  series  of  the  vertical  component 
of  velocity  on  geophone  #2  at  the  beginning  of  tape  RLAM-26  (see  Fig.  6-13).  Note 
the  significant  peaks  at  a  lag  of  about  0.8  sec  and  some  secondary  peaks  at  about 
0.35  sec  and  0.45  sec. 
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Figure  6-15:  Slip  function  for  part  of  the  time  series  shown  on  Fig.  6-13.  Note  the 
pattern  characteristic  of  stick-slip  motion,  when  the  process  of  slipping  stops  for  a 
time  (or  even  reverses)  and  then  continues  further  on. 
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Figure  6-16;  Spectrogram  of  vertical  velocity  on  geophone  #2  at  the  beginning  of 
the  tape  RLAM-26.  Note  that  the  regular  stick-slip  pattern  of  the  spectrum  changes 
after  about  12  sec. 
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When  two  surfaces  of  a  brittle  material  are  placed  together,  the  asperities 
on  the  surfaces  in  contact  become  locked  together.  If  the  normal  load  is 
high  enough  to  prevent  the  surfaces  from  lifting  up  over  the  irregularities, 
then  sliding  will  occur  when  the  locked  regions  fail  brittlely. 

The  evidence  for  the  type  of  motion  (stick-slip)  described  above  was  found  in 
the  geophone  data  from  the  ice  island  experimental  site.  On  Fig.  6-13  about  1.5  sec 
of  the  time  series  of  the  vertical  component  of  velocity  on  geophone  #2  are  shown. 
Note  several  pulses  similar  in  shape  at  approximately  1.1  sec,  1.4  sec  and  1.8  sec 
on  Fig.  6-13(a)  which  are  characteristic  for  a  stick-slip  type  of  motion.  Fig.  6-13(b) 
shows  the  details  of  the  pulse  at  1.4  sec  which  indicate  the  flexural  wave  arrival.  The 
autocorrelation  function  (see  Fig.  6-14)  confirms  the  similarity  between  the  pulses 
on  Fig.  6-13.  Fig.  6-15  depicts  the  slip  function  for  the  time  series  shown  on  Fig.  6- 
13.  Note  the  pattern  characteristic  of  stick-slip  motion,  when  the  process  of  slipping 
stops  for  some  time  (or  even  reverses)  and  then  continues  further  on.  On  Fig.  6-16  the 
spectrogram  of  vertical  velocity  on  one  geophone  is  shown;  observe  how  the  regular 
stick-slip  pattern  changes  after  about  12  sec. 

6.4  Event  parameter  estimation  from  the  combined 
results  of  spectral  processing  and  slip  history 
determination 

6.4.1  Event  parameter  definitions 

As  it  was  mentioned  above,  the  two  most  important  source  parameters  characterizing 
an  ice  event  are  seismic  moment  Mq  and  stress  drop  Act.  According  to  [22],  the 
seismic  moment  Mq  is  defined  as 


Mq  —  iiuA^ 


(6.5) 


170 


Figure  6-17:  Stress  drop  definition  according  to  [21].  Here  a\  is  the  initial  stress,  is 
the  final  stress,  is  the  effective  stress,  cr/  is  the  dynamic  frictional  stress,  and  a/r 
is  the  static  frictional  stress.  The  stress  drop.  Act,  refers  to  the  difference  between  a\ 
and  02- 

where: 

H  -  shear  strength  of  ice, 
u  -  average  slip, 

A  -  fault  area. 

To  determine  the  stress  drop  in  an  ice  event  we  need  to  consider  the  stress 
change  in  time  during  the  faulting  process  according  to  [21]  (see  Fig.  6-17).  In  order 
for  the  rupturing  to  begin  the  local  stress  must  rise  by  at  least  cTfrs  (the  difference 
between  the  frictional  stress  and  initial  stress).  A  slip  occurs  in  those  sections  of  fault 
where  this  condition  is  satisfied.  On  this  figure  the  stress  drop.  Act,  is  defined  as  the 
difference  between  the  initial  stress,  ai,  and  the  final  stress,  (72.  The  other,  often  used 
parameter,  effective  stress  denotes  the  difference  between  static  and  dynamic 
friction.  Sometimes,  instead  of  stress  drop  people  use  fractional  stress  drop  e  that  is 
defined  as  the  ratio  of  stress  drop  to  the  effective  stress  available  for  fracturing: 

Aa 

C7eff 
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According  to  [46]  the  relation  between  average  stress  drop  and  slip  is  given  by: 


=  Cfi  (1)  . 


(6.6) 


where 

{I  for  a  strike-slip  fault 

(6.7) 

for  a  dip-slip  fault 

In  Equation  (6.6)  u  is  again  the  average  slip,  and  L  is  the  characteristic  rupture 
dimension. 

6.4.2  Procedure  for  event  parameter  estimation 

In  Equations  (6.5)  and  (6.6)  all  parameters  on  the  right-hand  side  can  be  estimated 
using  the  results  from  spectral  and  slip  history  methods.  To  show  how  this  can  be 
done,  let  us  recall  from  the  previous  discussion  that,  according  to  [3],  the  corner 
frequency  is  inversely  related  to  the  dimensions  of  the  source  (or  source  duration  T): 


fo  =  I/ttT. 


Now,  to  obtain  the  length  of  a  crack  from  its  duration  T  we  need  to  make  an  assump¬ 
tion  regarding  the  rupture  velocity  Vr-  According  to  [31] 

Vr  «  0.63c„  (6.8) 

where  Cg  is  the  shear  wave  speed.  So  we  can  obtain  the  length  of  a  crack  L  using  the 
above  value  of  rupture  velocity: 


(6.9) 


A  value  of  p  —  920  kg/m^  for  the  density  of  sea  ice  is  used  for  calculations 
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based  on  data  provided  by  [38]: 

The  density  of  pure  ice  is  916.8  kg/m^.  However  the  density  of  sea-ice  may 
be  greater  than  this  last  figure  (if  brine  is  trapped  among  the  ice  crystals) 
or  less  (if  the  brine  has  escaped  and  gas  bubbles  are  present).  Values  from 
924  to  857  kg/m^  were  recorded  on  the  Norwegian  Maud  Expedition. 

The  values  of  other  important  constants  for  ice  are  given  below: 


Cg  =  1800  m/s 
//  =  3  X  10®  Pa 


So  finally,  using  Equation  (6.8)  and  (6.9)  and  values  of  the  constants  for  ice,  we  can 
estimate  the  event  parameters  from  the  following  relations: 


■gr  -A  -  Vyh 

Mq  =  ixuA  — 

TT/o 

Act  = 


(6.10) 

(6.11) 


where  h  is  the  ice  thickness  (2  m).  So,  we  see  that  Mq  and  Acr  can  be  estimated  using 
only  results  from  processing  by  both  spectral  and  slip  history  determination  methods 
and  values  of  appropriate  physical  constants  for  ice. 


6.4.3  The  results  of  estimation  of  event  parameters  and  com¬ 
parison  with  previous  estimates 

In  Table  6.1  results  of  parameter  estimation  using  both  spectral  and  slip  history 
determination  methods  and  Equations  (6.10)-(6.11)  are  given  for  several  events  from 
the  beginning  of  tape  RLAM-26.  In  this  table,  To  denotes  approximate  arrival  time 
for  the  event  in  seconds  from  the  beginning  of  tape  RLAM-26,  u  is  the  average  slip 
(after  correction  for  range  spreading),  /o  is  the  corner  frequency,  Mq  is  the  seismic 
moment,  and  Aa  is  the  average  stress  drop. 
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Event  # 

To,  sec 

w,  mm 

/o,  Hz 

Mo,  MNm 

Act,  kPa 

3 

37 

0.037 

7.8 

10.4 

1.54 

5 

39 

0.043 

8.8 

10.6 

2.00 

9 

58 

0.029 

8.4 

7.4 

1.27 

10 

65 

0.011 

7.5 

3.1 

0.42 

15 

145 

0.057 

17.6 

2.09 

16 

169 

0.056 

id 

15.3 

2.34 

Table  6.1:  Event  parameters  for  several  events  at  the  beginning  of  tape  RLAM-26. 
Note  that  units  of  seismic  moment  Mq  are  MNm  (i.e.  10®  N  x  m). 


These  results  give  the  following  range  for  stress  drop  values:  Aa  =  0.42  — 
—2.34  kPa.  The  range  of  values  for  seismic  moment  is  Mq  =  3.1 - 17.6  x  10®  Nm. 

In  the  thesis  of  C.  Stamoulis  [49],  stress  drop  ranges  from  110  Pa  at  /  =300  Hz 
to  2.25  X  10^  Pa  at  /  =  40  Hz  which  is  at  least  one  order  of  magnitude  higher  than  my 
upper  limit  of  values  for  stress  drop.  For  most  of  the  events  processed  by  C.  Stamoulis 
in  [49]  the  estimates  of  seismic  moment  Mq  lie  in  between  lO’^  and  10^®  Nm.  Only  very 
few  events  actually  go  below  10^  Nm  value,  but  nevertheless  some  overlap  between 
my  estimates  of  Mq  and  those  of  [49]  could  be  observed. 

To  compare  with  Chen’s  [4]  results  one  needs  to  be  aware  of  the  arithmetical 
error  on  page  105  of  [4]:  when  plugging  values  in  Equation  (4.42)  on  that  page  the 
actual  numerical  coefficient  obtained  is  9.9  x  10“^  not  5.9  x  10“^,  as  written  in  text. 
Hence  all  subsequently  derived  relations  should  be  multiplied  by  factor  1.68  to  correct 
for  that  error.  In  further  discussion,  first,  the  result  obtained  using  the  wrong  value 
will  be  shown,  then,  in  parenthesises,  the  result  using  the  correct  value  will  be  shown. 

The  stress  drop  range  is  predominantly  from  10  to  1000  Pa.  The  lowest 
frequency  data  point  is  at  approximately  70  Hz  with  a  value  of  stress  drop  being 
Act  =  5  kPa. 

After  adapting  notation,  Chen’s  regression  formula  for  the  dependence  of  stress 
drop  on  frequency  (Equation  (4.44)  on  page  105  of  [4])  takes  the  form 

Act  =  2.7  x  10“/"^"'^’’, 
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(the  correct  coefficient  should  be  4.6  instead  of  2.7). 

According  to  that  relation  for  frequency  /  =  10  Hz  we  have  Acr  =  2.7  kPa 
(4.6  kPa).  So  the  corrected  value,  4.6  kPa,  is  slightly  higher  than  the  upper  limit  of 
my  range  of  values  for  Aa  (2.34  kPa). 

To  estimate  the  seismic  moment  Chen  uses  an  octopole  strength  and  obtains 
the  following  expression  for  the  frequency  dependence  of  the  seismic  moment  (Equa¬ 
tion  (4.45)  on  page  106  of  [4]): 


Mo  =  4.1  X  lOV^-^ 

(the  correct  coefficient  should  be  6.9  instead  of  4.1).  Note  that  the  above-discussed 
error  propagates  into  the  results  for  the  seismic  moment. 

That  should  give  us  the  following  value  for  the  seismic  moment  at  /  =  10  Hz: 
Mq  =  13.0  X  10®  Nm  (actually,  21.8  x  10®  Nm).  Again,  the  corrected  value  is  slightly 
higher  than  the  upper  limit  of  my  range  of  values  of  the  seismic  moment  (17.6  x 
10®  Nm),  but  taking  into  account  the  variation  of  Chen’s  data,  one  can  say  that  her 
estimates  lie  within  the  range  of  event  parameters  estimated  in  this  thesis. 

Note,  that  the  ice  fracture  strength  lies  in  the  range  between  2  x  10®  and 
1  X  10®  Pa.  For  both  Chen’s  estimates  and  my  own,  the  difference  between  the  stress 
drop  in  ice  events  and  the  ice  fracture  strength  is  at  least  two  orders  of  magnitude. 
For  Stamoulis’  estimates  this  difference  is  at  least  one  order  of  magnitude.  This  is 
somewhat  similar  to  what  is  happening  in  seismology:  in  the  case  of  earthquakes, 
there  is  usually  a  difference  of  2-4  orders  of  magnitude  between  measured  stress  drops 
and  the  shear  strength  of  rock. 

6.5  Conclusions 

Using  spectral  methods  of  event  processing,  a  high  frequency  trend  rolling  off  as 
w”®  in  the  vertical  displacement  spectra  of  events  was  found  which  probably  could 
be  attributed  to  fault  nucleation.  In  some  of  the  events  at  even  higher  frequencies 
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a  second  trend,  was  observed.  Such  frequency  behavior  is  consistent  with  a 
commonly  accepted  model  in  earthquake  seismology  for  high  frequency  behavior  of 
the  displacement  spectrum  [21]. 

The  results  of  source  parameter  inversion  give  the  following  ranges  of  values: 

Stress  drop:  Aa  =  0.42 - 2.34  kPa. 

Seismic  moment:  Mq  =  3.1 - 17.6  x  10®  Nm. 

These  results  are  quite  close  to  estimates  by  Chen  for  the  marginal  ice  zone, 
but  somewhat  different  from  those  by  C.  Stamoulis  for  the  same  Arctic  region,  which 
probably  indicates  that  physical  conditions  on  the  “ice  island”  (where  all  the  data 
processed  in  this  chapter  were  collected)  are  closer  to  that  of  the  marginal  ice  zone 
than  to  those  usually  characteristic  of  the  Central  Arctic.  The  other  possibility  is 
that  the  kind  of  events  observed  on  the  “ice  island”  was  different  from  those  in  the 
vicinity  of  the  main  surveillance  array  which  was  the  area  covered  in  [49].  The  low 
level  of  pressure  on  the  lone  hydrophone  accompanying  the  RLAM  unit  during  the 
ice  events  discussed  above  seems  to  confirm  that  possibility. 

The  estimates  of  stress  drop  are  at  least  two  orders  of  magnitude  lower  than 
the  ice  fracture  strength  value  which  is  similar  to  earthquake  mechanics  case  (2-4 
orders  of  magnitude).  This  fact  probably  indicates  that  the  events  occured  on  the 
preexisting  fault  surfaces. 
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Chapter  7 


Conclusions  and  suggestions  for 
future  work 

7.1  Conclusions 

Using  the  Motion  Product  Detector  method  for  processing  multi-component  geophone 
data,  the  polarization  characteristics  of  elastic  waves  were  determined.  The  fractures 
processed  by  this  method  seemed  to  generate  vertically  polarized  shear  (SV)  waves 
for  the  most  part.  Their  concentration  in  the  vicinity  of  the  corner  of  the  new  ice 
ridge  and  the  character  of  the  particle  motion  in  the  corresponding  arrivals  indicate 
the  continuing  ridge  building  process  in  that  part  of  the  ice  island,  which  is  located 
approximately  4  km  East  of  the  base  camp. 

The  processing  of  data  from  an  ice  lead  approximately  2  km  North  of  the  base 
camp  provided  evidence  for  the  existence  of  a  new  phenomenon:  edge  waves,  waves 
propagating  along  a  newly  opened  lead.  The  waves  exhibit  a  quasi-periodic  behavior 
suggesting  some  kind  of  stick-slip  generation  mechanism  somewhere  along  the  length 
of  the  lead.  The  propagation  characteristics  of  these  waves  were  determined  using 
seismic  wavenumber  estimation  techniques.  In  the  low-frequency  limit  the  dispersion 
can  be  modeled  approximately  by  an  interaction  at  the  lead  edges  of  the  lowest  order, 
antisymmetric  modes  of  the  infinite  plate. 

The  qualitative  comparison  of  theoretical  radiation  patterns  with  radiation 
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characteristics  of  ice  events  suggests  a  strike-slip  fault  with  a  dip  angle  of  90°  as  the 
most  probable  source  mechanism  for  the  ice  events  processed  here,  which  is  somewhat 
to  the  contrary  of  J.S.  Kim’s  conclusion  [22].  On  the  basis  of  the  relative  smallness 
of  the  radiation  from  longitudinal  wave  compared  to  the  acoustic  mode  he  suggested 
that  dip-slip  with  dip  angle  6  0°  and  90°  and  strike-slip  with  dip  angle  5  ps  0°  are 

the  most  probable  source  mechanisms.  This  discrepancy  may  arise  from  the  fact  that 
most  of  the  events  processed  here  do  not  radiate  eflSciently  into  the  water. 

Using  spectral  methods,  the  high  frequency  trend  u}~^  in  the  vertical  dis¬ 
placement  spectra  of  events  was  found  which  probably  could  be  attributed  to  fault 
nucleation  processes.  In  some  of  the  events  at  even  higher  frequencies  a  second  trend, 
was  observed.  That  one  corresponds  to  the  high  frequency  behavior  of  the  dis¬ 
placement  spectrum  that  is  found  in  the  commonly  accepted  model  from  earthquake 
seismology. 

The  result  of  the  source  parameter  inversion  gives  the  following  ranges  of 

values: 

Stress  drop:  Aa  =  0.42 - 2.34  kPa. 

Seismic  moment:  Mq  =  3.1 - 17.6  x  10®  Nm. 

These  results  are  quite  close  to  estimates  by  Chen  for  the  marginal  ice  zone, 
but  somewhat  different  from  those  by  C.  Stamoulis  for  the  same  Arctic  region,  which 
probably  indicates  that  physical  conditions  on  the  ice  island  (where  most  of  the  data 
processed  in  this  thesis  were  collected)  were  closer  to  that  of  the  marginal  ice  zone 
than  to  those  usually  characteristic  of  the  Central  Arctic.  The  other  possibility  is  that 
the  kinds  of  events  observed  on  the  ice  island  were  different  from  those  in  the  vicinity 
of  main  surveillance  array  which  was  the  area  covered  in  [49].  This  difference  could 
be  due  to  the  active  ridge  building  processes  on  the  ice  island  which  were  absent  near 
base  camp  during  the  time  period  covered  by  the  data  processing  in  Stamoulis’  thesis 
[49].  The  low  level  of  pressure  on  the  lone  hydrophone  accompanying  the  REAM  unit 
during  the  above  discussed  ice  events  seems  to  confirm  that  possibility. 
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The  estimates  for  stress  drop  are  at  least  two  orders  of  magnitude  lower  than 
the  ice  fracture  strength  value  which  is  similar  to  the  case  of  earthquake  mechanics 
(2-4  orders  of  magnitude).  This  fact  probably  indicates  that  the  events  occured  on 
the  preexisting  fault  surfaces. 

7.2  Suggestions  for  future  work 

Additional  modeling  is  required  to  understand  more  completely  the  nature  of  the 
new  ‘edge  wave’  phenomenon  described  above,  in  particular  in  regard  to  achieving  a 
better  match  of  the  dispersion  behavior  with  the  model. 

The  other  direction  for  future  work  is  a  more  extensive  use  of  polarization 
methods  from  earthquake  seismology  for  geophone  data  processing. 


Appendix  A 


Inverse  for  an  over  const  rained 
system 


When  there  are  more  observations  than  unknown  model  parameters  (n  >  m)  the 
solution,  if  it  exists,  should  be  determined  by  a  least  squares  method.  Let  the  error 
for  each  data  set  be  e. 


e  =  Ax  —  y. 

The  squared  norm  is  obtained  from  a  product  of  and  e: 

=  e^e  =  (Ax  —  y)”^(Ax  —  y) 

(A.l) 

=  (x^A^  —  y"^)(Ax  —  y)  =  x^A^Ax  —  y^Ax  —  x^A^y  +  y^y 

To  minimize  the  squared  norm  we  differentiate  in  a  vector  space  with  respect  to  x  or 
x^  and  set  the  result  equal  to  zero: 

(A.2) 
(A.3) 


dx 

de^ 


—  x^A^A  —  y^  A  =  0 
=  A^Ax  —  A^y  =  0. 
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Both  equations  yield  the  same  result,  so  the  solution  for  the  estimate,  x*’,  could  be 
obtained  from  the  Equation  (A.3): 

[A'^A]  x**  =  A"^y  (A.4) 

x*^  =  ^  -A-^y-  (-A-S) 

It  is  necessary  that  the  matrix  [A^A]  be  non-singular: 

det(A'^A)  7^  0. 

The  matrix  [A'^A]  is  obviously  an  m  x  m  autocovariance  matrix,  while  A'^y 
is  a  cross  covariance. 
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